martes, 29 de noviembre de 2011

Modelos a la carta: El modelo gaussiano


Los modelos lineales generalizados (GLM) tienen tres componentes: (1) la estructura de errores; (2) el estimador lineal y; (3) la función de vínculo. La función de vínculo linealiza la relación entre la respuesta y el estimador lineal. La función de vínculo, junto con el estimador lineal, constituyen lo que lo algunos autores denominan el "modelo científico". La máxima versatilidad en la modelización se alcanza cuando uno ajusta modelos a la carta, fuera de las restricciones que los GLM imponen. Ello implica que debemos escribir nuestras propias funciones, para después ajustarlas a los datos. El proceso se completa en el contexto de la selección de modelos por criterios de información cuando comparamos toda una batería de modelos plausibles ajustados por métodos de máxima verosimilitud. En R esto puede hacerse por medio del paquete 'likelihood' escrito por Lora Murphy. El paquete todavía no está disponible en CRAN, aunque es de esperar que lo esté pronto. Mientras tanto, este paquete se puede descargar de la página del curso de "Modelos y métodos de máxima verosimilitud" de Charlie Canham. Este curso fue por cierto impartido en abril en Granada y he de decir que ha sido uno de los mejores cursos a los que he asistido... ¡gracias Charlie!

Bueno, a lo que iba... en lo que respecta a los modelos científicos. Podemos usar muchos tipos de modelos dependiendo de los datos que tengamos. No tienen por qué ser lineales. Una familia de modelos muy versátiles que permiten ajustar respuestas que aumentan primero y disminuyen después son los modelos gaussianos. La función gaussiana se expresa así:
y = a e ( ( X b ) 2 ( 2 c 2 ) ) y = a cdot e^(-(X-b)^2 over (2 cdot c^2))

Dicha función consta de 3 parámetros, que llamaremos 'a', 'b' y 'c', pero que podríamos haber denominado de cualquier otra forma. El parámetro 'a' indica el valor máximo que alcanza la curva. El parámetro 'b' determina la forma de la curva: monotónica, creciente, decreciente. El parámetro 'c' va a condicionar la pendiente de la curva (lo que en inglés se denomina sharpness).

Veamos ahora cómo cambios en los valores de estos parámetros determinan cambios en la forma de las curvas tipo. Tomemos un ejemplo en dónde vamos a modelar la abundancia de arácnidos en sistemas agrícolas (y) en función del tiempo (X). La idea con el siguiente código en R es que nos familiaricemos con esta función y sus parámetros.

Si modificamos el parámetro 'a' podemos obtener algunas las siguientes curvas:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=125, y=1025, "a=1000")
a <- 800
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=125, y=825, "a=800")
a <- 600
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=3)
text(x=125, y=625, "a=600")
a <- 400
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=125, y=425, "a=400")
a <- 200
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=125, y=225, "a=200")
Modificando el parámetro 'b' obtenemos cambios en la forma de la curva gaussiana como se observa en la siguiente gráfica:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=125, y=1025, "b=125")
b <- 200
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=200, y=1025, "b=200")
b <- 275
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=240, y=825, "b=275")
b <- 50
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=50, y=1025, "b=50")
b <- -25
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=15, y=825, "b=-25")
Cambios en el parámetro 'c' van a provocar cambios en las pendientes de la curva, como se observa a continuación:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=45, y=825, "c=0.01")
c <- 0.001
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=125, y=1025, "c=0.001")
c <- 0.25
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=125, y=150, "c=.025")
c <- 0.1
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=90, y=250, "c=0.1")
c <- 0.02
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=40, y=400, "c=0.02")

4 comentarios:

cjgb dijo...

¿Has probado las funciones mle de stats4 o fitdistr de MASS? ¿Te valen para lo que intentas? ¿O me pierdo algo?

Luis Cayuela dijo...

Realmente no intento nada... era un caso para ilustrar lo qué significan los parámetros de la función gaussiana. Realmente cuando utilizo este tipo de modelos, suelo complicarlos con covariables y otros términos que se añaden a la función y ajusto los parámetros por métodos de MLE con el paquete 'likelihood'. Pero no era este el caso...

Gracias en cualquier caso por el comentario y por las sugerencias.

cjgb dijo...

Es que como había visto dos entradas tuyas sobre asuntos "pre-MLE" (el de la obtención del MLE por fuerza bruta y éste) pensé que tal vez sería un colofón mencionar esas funciones para obtener el MLE numéricamente, aunque fuese como una nota.

Al fin y al cabo, no son funciones que mucha gente conozca.

David dijo...

Todo lo relacionado con la docencia me importa mucho ya que me dedico a la educación y por eso trato de ayudar a los niños. Muchas veces fui a fravega para conseguir lo que necesito sobre tecnología y por eso trato de utilizarla vinculada a la educación

Buscar entradas