FUNCIONES
Creación de funciones
Cómo vimos, una función es una serie de sentencias organizadas en
bloque para realizar una tarea específica. R tiene muchísimas funciones
y los paquetes externos son basicamente colecciones de funciones
organizadas en una temática o tipo de análisis. Además de estas
funciones pre-establecidas (sum(), max(),
min(), mean()…), podemos crear nuestras
propias funciones.
Las funciones se crean usando la función function() (que
fue creada por generación espontanea). La sintáxis básica sería:
function_name <- function(arg_1, arg_2, ...) { Function body (tareas que queremos realzar) }
Las funciones tienen diferentes partes o componentes:
Nombre de la función: es la manera en la que vamos a llamarla para usarla. Al no ser una función pre-establecida tenemos que habrela cargado (ejecutado el código) para que esté en el entorno de trabajo.
Argumentos: Los argumentos son marcadores de posición. Cuando llamamos a una función le damos esos valores para poder realizar la tarea exactamente como la queremos realizar. Los argumentos son opcionales (pueden no existir) y pueden tener valores por defecto.
Cuerpo de la función: Contiene el bloque con la colección de sentencias que define qué hace la función.
Salida: Son los resultados de la función. Toda función debe tener al menos un resultado. Si una función entrega varios tipos de objetos se acostumbra a organizarlos en una lista que puede manejar los diferentes tipos de objetos.
Veamos un ejemplo
suma_dos<-function(n1,n2){
n1+n2
}
suma_dos(3,4)## [1] 7Si queremos ser prolijos (y seguir la lógica), esto en verdad debería haberse escrito:
suma_dos<-function(n1,n2){
salida<-n1+n2
return(salida)
}
suma_dos(3,4)## [1] 7Pero bueno, usemos la versión abreviada por comodidad…
Los argumentos pueden tener un valor por defecto, por ejemplo
calculo<-function(a,b,c=1) {
a*b^c
}
calculo(a=2, b=3)## [1] 6calculo(a=2, b=3, c=3)## [1] 54Cuando hay mas de un argumento se pueden proveer los valores sin nombrarlos (toma el orden en el que aparecen) o podemos nombrarlos y ya no importa el orden, por ejemplo
calculo2<-function(a,b,c) {
a*b^c
}
calculo2(2,3,3)## [1] 54calculo2(c=2, b=3, a=3)## [1] 27Si un argumento no está (y no tiene valor por defecto), entonces la función va a dar error
calculo2(a=4, c=5)## Error in calculo2(a = 4, c = 5): el argumento "b" está ausente, sin valor por omisiónA veces queremos que la salida sea mas de un resultado. Para lograrlo tenemos que especificar qué queremos que sea la salida, ya que si no mostrará solamente la última tarea realizada.
calculo3<-function(a,b,c) {
a*b
a+b
b^2
c+3
}
calculo3(a=4, b=5, c=6)## [1] 9calculo3<-function(a,b,c) {
print(a*b)
print(a+b)
print(b^2)
print(c+3)
}
calculo3(a=4, b=5, c=6)## [1] 20
## [1] 9
## [1] 25
## [1] 9#otra posibilidad usando return()
calculo3<-function(a,b,c) {
s1<-a*b
s2<-a+b
s3<-b^2
s4<-c+3
return(c(s1,s2,s3,s4))
}
calculo3(a=4, b=5, c=6)## [1] 20 9 25 9Sin embargo la evaluación de la existencia de valores para los argumentos es “lazy” y solo se realiza al momento de usar ese argumento:
calculo3<-function(a,b,c) {
print(a*b)
print(a+b)
print(b^2)
print(c+3)
}
calculo3(a=4, b=5)## [1] 20
## [1] 9
## [1] 25## Error in calculo3(a = 4, b = 5): el argumento "c" está ausente, sin valor por omisiónPodemos complejizar un poco y hacer la salida mas “amigable” para el
usuario. Supongamos que queremos calcular el producto entre dos valores
y que la salida sea un texto (usamos la función cat() para
unir texto con los valores que queremos que muestre):
producto <- function(a, b) {
produ <- a * b
cat("El producto entre", a, "y", b,
"es", produ)
}
producto(a=3, b=9)## El producto entre 3 y 9 es 27INTEGRANDO FUNCIONES Y LOOPS
Veamos un ejemplo integrando lo que vimos sobre loops y funciones
Usemos datos de tallas de Sander vitreus, un pez dulceacuícola que se distribuye en los lagos del norte de Estados Unidos y de Canada. En el archivo peces.txt se presentan datos de la edad (años) y la longitud (mm) de peces hembra capturados en 3 lagos dentro del rango de distribucion de la especie. Podemos usar la funcion de crecimiento de von Bertalanffy para modelar el crecimiento de esta especie a partir de los datos. \[L_{Edad} = L_{∞} ∗ (1 − exp(−k ∗ (Edad − t_0)))\] En este modelo \(L_∞\) representa la longitud asintotica, k la tasa de crecimiento somatico, y \(t_0\) la talla esperada a la edad 0 (en general un valor negativo). Tenemos entonces una formula teórica de cómo debería ser el aumento de talla a lo largo del tiempo, y tenemos los datos empíricos. Empecemos cargando los datos y mirando su estructura
crecimiento<-read.table("peces.txt", header=T)
head(crecimiento)## Edad Longitud
## 1 0 188
## 2 0 190
## 3 0 191
## 4 0 196
## 5 0 173
## 6 0 176Veamos como se ve la relación entre edad y tamaño
plot(crecimiento)
Podemos definir una función que aplique von Bertalanffy para calcular el largo de un pez dada su edad y ciertos valores de parámentros. Vamos a usar x para la edad y par para los parametros (osea par debería contener 3 valores, uno para \(L_∞\), otro para k, y otro para \(t_0\)).
vonberta<-function(x, par)
{
par[1]*(1-exp(-par[2]*(x[, 1]-par[3])))
}Como no sabemos los valores de los parámetros, podríamos aproximarlos a ojo. El largo máximo es medianamente fácil porque es la asintota. Pongamos un valor un poco mayor que el largo maximo observado y los demas parámetros los inventamos.
max(crecimiento$Longitud)## [1] 739vb<-vonberta(x=crecimiento, par=c(900,0.1,0))vb tiene, para cada observación que teníamos en la planilla de datos, los valores esperados de longitud según von Bertalanffy para ese individuo.
vb[1:30]## [1] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [9] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## [17] 85.64632 85.64632 85.64632 85.64632 85.64632 85.64632 85.64632 85.64632
## [25] 85.64632 85.64632 85.64632 85.64632 85.64632 85.64632Veamos como nos fue eligiendo a ojo los valores de parámetros. Primero graficamos los datos y le superponemos la curva calculada con la función
plot(crecimiento)
lines(crecimiento[,1], vb)
Podríamos cambiar los valores de parámetros y repetir a ver si nos va mejor.Intentemos achicando la talla al nacer.
vb2<-vonberta(x=crecimiento, par=c(900,0.1,-2))
plot(crecimiento)
lines(crecimiento[,1], vb2)
Nos fue mejor. Pero tendríamos que encontrar una forma de definir “mejor” mas objetiva que el ojo. Una manera intuitiva de definir “mejor” es tener en cuenta cuan lejos quedan los valores reales (los valores observados) con respecto a los valores estimados por la curva. Esto es; la famosa suma de cuadrados del error. Ya que estamos definiendo funciones, podemos definir una función que justamente compute la suma de cuadrados del error.
ajuste_vb<-function(par){
yesp=vonberta(x=crecimiento, par=par)
y=crecimiento[,2]
sum((y-yesp)^2)
}Veamos el ajuste de los dos intentos que hicimos:
ajuste_vb(par=c(900, 0.1, 0))## [1] 17528726ajuste_vb(par=c(900, 0.1, -2))## [1] 1473180Podríamos pasarnos horas intentando combinaciones de valores de
parámetros diferentes e ir viendo el ajuste para quedarnos con los
parámetros que lo maximicen (minimicen la suma de cuadrados del error).
Existe una función, optim(), que realiza las iteraciones
por nosotros y calcula los valores óptimos de parámetros. Esta función
necesita un valor inicial de parámetros. Usemos los que vimos que
parecían mejores.
par_ini=c(900, 0.1, -2)
AjuVB=optim(ajuste_vb, par=par_ini, method="BFGS")
str(AjuVB)## List of 5
## $ par : num [1:3] 729.895 0.166 -1.671
## $ value : num 621261
## $ counts : Named int [1:2] 132 24
## ..- attr(*, "names")= chr [1:2] "function" "gradient"
## $ convergence: int 0
## $ message : NULLAjuVB$par## [1] 729.8945522 0.1660436 -1.6706616Tenemos los parámetros del mejor ajuste, podemos usarlos para graficar.
VB_posta<-vonberta(x=crecimiento, par=AjuVB$par)
plot(crecimiento)
lines(crecimiento[,1], VB_posta)
Como los datos provienen de un muestreo, teoricamente existen valores de estos parámetros que describen el crecimiento de estos organismos (los parámetros poblacionales) y nosotros, mediante el muestreo, estamos tratando de estimarlos. Sabemos que nuestras estimaciones van a tener cierto nivel de incertidumbre (dado por usar un subconjunto de la población). En estadística paramétrica poddríamos usar la teoría para calcular esta incertidumbre. En este caso, que no contamos con teoría podríamos usar técnicas de remuestreo para calcularla. La idea es, a partir de nuestros datos, generar diferentes muestras, que difieran de la original, para ver cuanto difieren los valores estimados de los parámetros. Si la insertidumbre es poca, entonces la variabilidad en la estimación de los parámetros va a ser baja. Las diferentes muestras las podemos obtener remuestreando los datos con reposición. Primero nos fijamos cuantos datos tenemos y creamos un indice que sea el remuestreo de esos datos. Como es un muestreo con reposición, terminamos con el mísmo número de datos pero algunos pueden repetirse y otros faltar. Usamos ese indice para crear nuestros datos remuetreados.
nrow(crecimiento)## [1] 707indi=sample(1:707,rep=T)
sort(indi)[1:20]## [1] 3 4 6 6 7 7 7 8 11 11 12 12 12 12 13 13 13 15 16 18remuestreo<-crecimiento[indi,]
head(remuestreo)## Edad Longitud
## 568 5 551
## 246 2 350
## 462 4 455
## 567 5 544
## 359 3 392
## 530 4 484Una vez que tenemos la lógica de cómo crear las muestras, podemos usar un loop para realizar el procedimiento n veces y, en cada ciclo del loop, quedarnos con los valores estimados de los parámetros.
salida=NULL
for(i in 1:100){
indi=sample(1:707,rep=T)
l_ajuste_vb<-function(par){
yesp=vonberta(x=crecimiento[indi,], par=par)
y=crecimiento[indi,2]
sum((y-yesp)^2)
}
l_AjuVB=optim(l_ajuste_vb, par=par_ini, method="BFGS")
salida<-rbind(salida, l_AjuVB$par)
}
head(salida)## [,1] [,2] [,3]
## [1,] 740.8841 0.1605284 -1.691123
## [2,] 708.4244 0.1753676 -1.637223
## [3,] 758.5459 0.1522835 -1.800175
## [4,] 723.3834 0.1733373 -1.569389
## [5,] 767.1207 0.1479693 -1.829733
## [6,] 707.0861 0.1764174 -1.614337Ahora que tenemos los valores estimados podemos graficarlos, o calcular intervalos de confianza. Veamos, por ejemplo, la distribución de las estimaciones de la longitud asintotica.
hist(salida[,1])
c(mean(salida[,1]), quantile(salida[,1], 0.025), quantile(salida[,1], 0.975))## 2.5% 97.5%
## 732.4098 705.4068 768.4150