ANÁLISIS DE DESTINOS CONOCIDOS CON COVARIABLES INDIVIDUALES en RMARK
Algunas de las aplicaciones más interesantes e importantes de los modelos estadísticos van más allá de la estimación de parámetros básicos como la supervivencia, y profundizar en describir y explicar la heterogeneidad en la supervivencia y otras tasas vitales (idealmente en situaciones experimentales, como cuando los factores que han sido manipulados y los sujetos han asignados a grupos de tratamiento). Programas MARK y el paquete RMark tienen enorme capacidad para la inclusión de este tipo de modelos, en el supuesto de que la toma de muestras o el diseño experimental apropiado se ha seguido.
Ilustro esta situación con un ejemplo de conjunto de datos de telemetría de los Patos Negros Americanos (Conroy et al. 1989. Journal of Wildlife Management 53:99-109). Los datos están contenidos en una trama de datos incorporada en RMark, que podemos acceder a través de los comandos
> require(RMark)
> data(Blackduck)
Se lista las primeras líneas de datos que tenemos ...
> head(Blackduck )
ch BirdAge Peso Wing_Len condix
ch BirdAge Weight Wing_Len condix
1 1100000000000000 1 1.16 27.7 4.19
2 1011000000000000 0 1.16 26.4 4.39
3 1011000000000000 1 1.08 26.7 4.04
4 1010000000000000 0 1.12 26.2 4.27
5 1010000000000000 1 1.14 27.7 4.11
Cada línea es una historia de captura (en este caso para los 5 primeros animales) en formato "LDLD", en la que "L" es sinónimo de observación en un intervalo de muestreo (1 = observado, 0 = no observado) y "D" significa supervivencia (0) o la muerte (1) en un intervalo. Esto es seguido por 4 covariables individuales: edad (1 = adultos, 0 = subadulto), masa corporal en kg, cuerda del ala en cm, y un índice de condición = masa / ala * 100. Obsérvese que en la trama de datos de BirdAge RMark se trata como una variable de factor (de modo no numérico), que puede producir resultados impredecibles en RMark. La expresión
>Blackduck$BirdAge=as.numeric(Blackduck$BirdAge)-1
la convierte a una covariable numérico (valores de 0 o 1).
Los datos se procesan y se asocian a la estructura del modelo destino conocido por
> bduck.processed=process.data(Blackduck,model="Known")
y una matriz de diseño general se construye como
> bduck.ddl=make.design.data(bduck.processed)
Finalmente, especifico una covariable tiempo específico, que es el número de días en una semana (ocasión) en el que la temperatura mínima fueran <0 C (congelación); esto se incorpora en la matriz de diseño.
> bduck.ddl$S$min=c(4,6,7,7,7,6,5,5)
Ahora podemos proceder a desarrollar varios modelos de interés. La sintaxis general para hacer esto es bastante simple, ya que sólo hay un único tipo de parámetro (S = supervivencia) que estamos modelando. Por ejemplo, el modelo que implica la constante S puede ser creado por los estados
>S.dot=list(formula=~1)
>mod1=mark(data=bduck.processed,ddl=bduck.ddl,model.parameters=list(S=S.dot))
Un modelo específico de tiempo se puede construir (con el "tiempo" indicador incorporado de parámetro de la matriz de diseño) como
>S.time=list(formula=~time)
>mod2=mark(data=bduck.processed,ddl=bduck.ddl,model.parameters=list(S=S.time))
Una serie de modelos se puede crear y ejecutar mediante la función mark.wrapper. Aquí está el código que se ejecuta estos modelos 2 y varias más, y genera una tabla resumen de los valores de la AIC. Los resultados parecen favorecer el modelo con la temperatura (min) como covariable, con el coeficiente que tiene un valor de -0,61, lo que indica relación negativa entre los días de congelación y la supervivencia. El código también produce un archivo de entrada que se puede ejecutar en la interfaz gráfica de usuario MARK.
La función de RMark model.average() se puede utilizar con el objeto bduck.result para calcular estimados de los modelos de promedio poderados de los "parámetros reales" (supervivencia semanal para semanas 1-8).
> model.average(bduck.results)
Tenga en cuenta que cuando las covariables individuales están involucrados (ya que son por lo menos con 2 modelos) los estimados reales de los parámetros realemente son predicciones en el valor de covarianza media, mientras que covariable específica a tiempo asumir el valor ocasión específica para la covariable. Por ejemplo, para el modelo S (min + BirdAge * Peso) que predeciría S (3) para los adultos mediante el establecimiento BirdAge = 1, min = 7, y con el peso promedio de los adultos. Predicción más flexible es habilitar mediante una función () covariate.predictions, tal como se describe a continuación.
Predicción mediante modelos de covarianza
A menudo, una vez que tenemos un estimado de los parámetros de un modelo de covarianza estaremos interesados en hacer predicciones sobre la respuesta a varios niveles de la covariable, incluidos los que pueden no haber sido observado aún. Si tenemos un modelo logit lineal de la forma
logit (S) = intercepto + pendiente * covariable
y tenemos valores estimados de la intersección y la pendiente (de nuestro análisis), podríamos simplemente conecte los diversos niveles de la covarianza y producir preditions. Sin embargo, esto se hace un poco engorroso para un par de razones: Normalmente le queremos valores estimados de la varianza y tal vez los intervalos de confianza en nuestras predicciones, y estos son computacionalmente más difícil que sólo los cálculos "plug-in" .
Además, la mayoría de las veces tendremos valores estimados bajo varios modelos, que luego producen correspondientemente diferentes predicciones dadas las valores de covarianza (con algún supuesto predecir ningún efecto de la covariable).
Afortunadamente, Jeff Laake ha escrito una muy buena función covariate.predictions () que se ocupa de estas cuestiones de una manera muy general. La función lee una trama de datos de los resultados (como los bduck.results anteriores), que a su vez todas las tiendas de los valores estimados de los parámetros del modelo en cada modelo, los valores de AIC, etc. Lee una trama de datos construida por el usuario que contiene los valores de las covariables que queremos predecir con. Este puede ser un rango de una sola covariable, o combinación de múltiples covariables.
La función crea un objeto con las predicciones del modelo promediado y la estructura de varianza-covarianza en cada valor de covarianza (o combinación de valores).
La ilustro con el ejemplo de Patos Negros. Supongamos que ya hemos creado el objeto de bduck.results del anterior. A continuación, he seleccionado una gama de valores de masa sobre la que predecir. Como notaron estos podrían ser cualquier cosa, por simplicidad estos son espaciados uniformemente en todo el rango de datos.
> ########################covariate predictions
> #pick a range of masses to predict over. Just use data range for now
> minmass=min(Blackduck$Weight)
> maxmass=max(Blackduck$Weight)
> mass.values=minmass+(0:30)*(maxmass-minmass)/3
Las predicciones del modelo promediada se consiguen por un solo muy simple llamada a la función covariate.predictions ()
> #compute model-averaged predictions of S at each mass value
> Sbymass=covariate.predictions(bduck.results,data=data.frame(Weight=mass.values),indices=c(1))
>
Si queremos, podemos trazar la predicción y el intervalo de confianza de todo el rango de valores de covarianza.
> minmass=min(Blackduck$Weight)
> maxmass=max(Blackduck$Weight)
> mass.values=minmass+(0:30)*(maxmass-minmass)/30
>
> #compute model-averaged predictions of S at each mass value
> Sbymass=covariate.predictions(bduck.results,data=data.frame(Weight=mass.values),indices=c(1))
> #plots
> plot(Sbymass$estimates$covdata, Sbymass$estimates$estimate,
+ type="l",lwd=2,xlab="Mass(kg)",ylab="Survival",ylim=c(0.8,1))
> lines(Sbymass$estimates$covdata, Sbymass$estimates$lcl,lty=2)
> lines(Sbymass$estimates$covdata, Sbymass$estimates$ucl,lty=2)
También pudimos ver cómo la predicción se ve en una superficie más extrema (500 g patos - probablemente muertos a los patos 2kg muy grasa!)
> #a more extreme prediction
> minmass=0.5
> maxmass=2.0
> mass.values=minmass+(0:30)*(maxmass-minmass)/30
>
> #compute model-averaged predictions of S at each mass value
> Sbymass=covariate.predictions(bduck.results,data=data.frame(Weight=mass.values),indices=c(1))
> #plots
> plot(Sbymass$estimates$covdata, Sbymass$estimates$estimate,
+ type="l",lwd=2,xlab="Mass(kg)",ylab="Survival",ylim=c(0.8,1))
> lines(Sbymass$estimates$covdata, Sbymass$estimates$lcl,lty=2)
> lines(Sbymass$estimates$covdata, Sbymass$estimates$ucl,lty=2)
>
Usted puede ver por sí mismo que, si bien las predicciones a 500 g y 940g (extremo inferior del rango de datos) no son enormemente diferentes, el intervalo de confianza inferior para las predicciones fuera del rango de datos se vuelve muy amplia (de la carta). Esto no debería ser una gran sorpresa y simplemente refuerza la advertencia acerca de la "predicción fuera del diseño."
Comentarios adicionales
Predicción de covariables
En lo anterior, la nota que indice = c (1) se refiere a los índices de los parámetros en los "todos diferentes" índices de parámetros (el índice de los parámetros reales), que para esta estructura de datos son 1-8 para los 8 ocasiones el tiempo. Los modelos que no incluyen un efecto de tiempo harán las mismas predicciones para todos los índices de los parámetros, porque es necesaria usar sóla el primero. En este ejemplo estamos promediando sobre algunos modelos que incluyen efectos en tiempo, por lo que la predicción promedio de modelo anterior referiremos sólo a la primera ocasión; si se necesitan predicciones para todas las ocasiones, esos índices tendrían que ser incluidos (índices = 1:8). Sin embargo, dado que el efecto covariable es el mismo en todas las ocasiones en los modelos donde hay un efecto, no hay realmente ninguna razón para hacer esto para todos los índices, especialmente si el interés se centra en el efecto comparativo del índice. Vamos a discutir las matrices de índice de parámetro (PIMS) mucho más a fondo, cuando lleguemos a los modelos CJS.
Puede comparar los resultados de la covariate.prediction () para enchufar los valores en el modelo de covarianza. Tome el caso simple de la edad + modelo de peso. En primer lugar os muestro el código que se saca de las estimaciones a partir del modelo, se conecta en valores de covarianza, y produce predicciones. Entonces hago lo mismo usando la función () covariate.prediction, ilustrando que la entrada para esa función puede ser un modelo único, - donde tenemos una predicción bajo sólo ese modelo - o una lista de los modelos (resultados) - donde tenemos un modelo de predicción promedio.
> #bird age +weight model
> betas<-bduck.results$S.BirdAge.Weight$results$beta$estimate
> #covariate inputs for intercept =1 age=1 and weight =1.2
> x<-c(1,1,1.2)
> #compute prediction
> SpredA<-expit(sum(betas*x))
> SpredA
[1] 0.9289718
>
> #using the covariate prediction function & just one model
> SpredB<-covariate.predictions(bduck.results$S.BirdAge.Weight,data=data.frame(BirdAge=1,Weight=1.2),indices=c(1))
> SpredB$estimate$estimate
[1] 0.9289718
>
Parámetro derivado
Los parámetros reales de supervivencia S modelados en RMark refieren a las probabilidades de sobrevivir cada intervalo del estudio, en el ejemplo de Patos Negros cada uno de los 8 períodos semanales. Aunque no forma parte explícitamente de la modelo de destino conocido, el interés se centra a menudo en la supervivencia al final del estudio, que por definición es el producto de S (j) sobre los j = 1, ..., n períodos, de modo S (1) * S (2) * S ....... (8) si la supervivencia es variable en el tiempo o S ^ 8 si no lo es. Esta cantidad también está disponible como un parámetro derivado en cada modelo, por ejemplo
> bduck.results$S.dot$results$real$estimate^8
[1] 0.5922517
muestra cómo se puede calcular a partir del valor del parámetro real y
> bduck.results$S.dot$results$derived
estimate se lcl ucl
1 0.5922518 0.07313623 0.4451435 0.7244965
proporciona este valor directamente a lo largo con su IC se y.
Datos de éxito de los nidos
La estructura de formato de datos y el modelo anterior es adecuada para el tipo usual de los datos de "destinos conocidos" encontrados en los estudios de radiotelemetría. Estudios de éxito de los nidos de las aves incluyen una estructura de datos un poco diferente, relacionado con la forma en que se determinan los destinos de nidos. Los datos de cada nido implican (1) el día en que el nido fue encontrado por primera vez, (2) último día presente = último día en que se observaron la presencia en el nido de polluelos, (3) el último día en que el nido fue detectado, y (4) una indicador de destino nido 0 = tramó 1 = depredado. También se especifica una columna de frecuencia; lo cual suele ser 1. Además, covariables (-nido específico) individuales pueden ser grabadas.
He incluido guiónes de R para ejecutar análisis 2 ejemplos de éxito de nidos. El primero es un simple análisis de los datos sobre 18 nidos de tero norteamericano, que se ejecuta una serie de modelos de ajuste términos de efectos en tiempo polinomial (así modelar la variación en la tasa de supervivencia diaria durante la temporada). El segundo es un estudio grande (565 nidos) del éxito de nidos de patos silvestres, en la que el éxito se modela como una función de varias covariables individuales (-nido específico) y solares se forman mostrando predicciones entre la tasa de supervivencia diaria y la edad del nido, y DSR y el hábitat características. Ver RMark ayuda para obtener más documentación sobre estos ejemplos.
Siguiente - Ejercicios