Bases del modelos del MARK - las matrices de PIM y Diseño
Independientemente de si usamos RMark o MARK, tenemos que estar familiarizado con algunos conceptos básicos para la estructuración de un modelo. Hay 2 elementos esenciales: la matriz de índices de los parámetros (PIM) y la matriz del diseño (DM).
PIM
La PIM es básicamente una lista de parámetros de diferentes tipos, con un número de indexación asociado a cada uno. Los tipos de parámetros y sus índices son determinados por (1) el tipo básico de los datos (por ejemplo, CMR para poblaciones cerradas vs. abiertas ), y (2) las dimensiones de los datos (el número de ocasiones de muestreo y los grupos o estratos). Vamos a ilustrar para la estructura básica del modelo para CMR para una población cerrada y un solo grupo (estrato) con 5 ocasiones de captura. Este modelo y la estructura de datos tendría un PIM con 3 tipos de parámetros diferentes : las probabilidades de captura (p), las probabilidades recaptura (c) y la abundancia (N). Una PIM básica que permite este tipo de parámetro 3 para tener valores distintos unos de otros y más ocasiones sería numerar los parámetros:
01:05 P1-P5
C1-C-4 06:09
N 10
Es decir, son 5 probabilidades de captura distintas permitidos en los 5 ocasiones, seguido de 4 probabilidades de recaptura (tenga en cuenta que mientras estos se clasifican de c1-c4 que se refieren a la probabilidad de la recaptura de un animal marcado previamente capturado en ocasiones 2-5). Finalmente, el último de estos parámetros es la abundancia (N), que ya estamos asumiendo que es cerrada, obviamenta no varía. Tenga en cuenta que es importante que los parámetros tienen índices distintos, es decir normalmente no vamos a permitir que los parámetros de los diferentes tipos de "compartir" un índice. Una excepción a esto se verá más adelante, donde permitimos p y c para compartir índices (ya que en algunos modelos queremos que estos parámetros sean los mismos).
Matriz de diseño (DM)
El DM opera en una estructura modelo básico definir por el PIM, y crea modelos particulares en supuestos específicos. por ejemplo
especifica un modelo en el que p y c son constantes en el tiempo pero diferentes entre sí (y, por supuesto, a partir de N); este modelo tiene 3 parámetros. En comparación, el modelo
indica que la captura y recaptura son las mismas (p = c) y constante, y tiene 2 parámetros. Tenga en cuenta que las columnas de la matriz de diseño terminan correspondiente al número de parametros que estamos tratando de estimar en el modelo.
En la práctica vamos a crear modelos usando el lenguaje de script en RMark, como se verá en los ejemplos siguientes. RMark combina esencialmente los pasos PIM y DM en un solo paso, eliminando los índices de los parámetros si estos no se utilizan en última instancia, en el diseño. Para ilustrar, si construimos un modelo con los efectos del tiempo y de comportamiento, podríamos utilizar un script similar a este para definir la probabilidad de captura
(time + c)
y el modelo tendría que seguir pasos de los índices de tiempo de 1-5 para la captura y recaptura. Si por el contrario sólo estábamos interesados en un efecto de la conducta, un modelo como
(~ 1 + c)
se utilizaría, y no habría necesidad de realizar un seguimiento de los índices de tiempo (ya que el efecto no implican tiempo). Esto será más claro con ejemplos específicos, más adelante. El ~ 1 indica que estamos reemplazando los efectos del tiempo para ocasiones 1-5 con una intercepción.
También tenga en cuenta que podríamos, en teoría, crear un modelo en el que el número de parámetros estimados igualó la lista original en el PIM (en el ejemplo anterior 10), pero, de hecho, que casi nunca va a ser el mejor modelo, y por lo general va a ser mejor para la construcción de una serie de modelos alternativos más simples.
Configuración de un análisis en RMark
El primer paso en un análisis RMark (o MARK) está recibiendo los datos en el formato correcto, que para RMark estará en una trama de datos especial. Por el momento vamos asumir que tenemos nuestros datos formateados, y utilizar un ejemplo en lata para ilustrar la construcción de modelos y estimación. Luego regresaremos y mostramos cómo poner los datos en RMark.
Primero tenemos que cargar la biblioteca RMark a buscar un ejemplo de datos. Los siguientes comandos traer un ejemplo clásico de Edwards y Eberhardt.
> library (RMark)
> data (edwards.eberhardt)
El comando
> edwards.eberhardt
muestra los datos. Usted verá que los datos están en una matriz de 76 filas y 18 columnas, que contiene unos y ceros. Las filas son las historias de captura de animales individuales (76 en total) y el 1/0 indica captura (1) o no captura (0) en cada una de las 18 ocasiones.
La sintaxis básica para leer los datos, creer un modelo, y obtener valores estimados es bastante simple, y se ve de forma genérica como éste
model_object <-mark (data, model, model.parameters)
"datos" nos dice dónde están los datos, "model" especifica la estructura del modelo (en este caso, cerrado CMR), y model.parameters indicará el modelo específico que queremos construir. El comando
mb <-mark (edwards.eberhardt, model = "Closed")
en realidad crear un modelo simple (pero no el más simple) para esta estructura y guardar el modelo y cualquier resultado en el objeto "mb". Por defecto, el modelo básico tiene un sólo parámetro para cada tipo de parámetro (p, c, y N), así que no hay variación en la captura a través del tiempo. Hay por lo tanto 3 parámetros. Podemos ejecutar este y ver algunos resultados de inmediato:
mb<-mark(edwards.eberhardt,model="Closed")
.........................
...............
> mb$results$real
estimate se lcl ucl fixed note
p g1 t1 0.0868518 0.0215415 0.0528969 0.1393948
c g1 t2 0.0809816 0.0095560 0.0641189 0.1017966
N g1 93.8918520 10.8161480 81.9921320 129.4231100
> mb$results$AICc
[1] 381.5498
Tenga en cuenta por la forma en que la formulación explícita de este modelo da exactamente los mismos resultados
pdot.cdotshared=list(formula=~1+c,share=TRUE)
mb_e<-mark(edwards.eberhardt,model="Closed", model.parameters=list(p=pdot.cdotshared))
Así, bajo este modelo, obtenemos valores estimados de N de aproximadamente 94, de la captura de p = 0,087, y de la recaptura de 0.081, por lo que podemos concluir que no existe (alta) evidencia de trampa timidez.
Además, un vez se crea un modelo podemos volver al su objecto y confirmar (por ejemplo) los PIMS el la matriz del diseño. Por ejemplo,
> PIMS(mb,"p")
group = Group 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> PIMS(mb,"c")
group = Group 1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
> PIMS(mb,"N")
group = Group 1
[,1]
1 3
> mb$design.matrix
p:(Intercept) c:(Intercept) N:(Intercept)
p g1 t1 "1" "0" "0"
c g1 t2 "0" "1" "0"
N g1 "0" "0" "1"
>
Volveremos a este punto más tarde por modelos más complejos.
Pero, ¿realmente necesitamos aún esta cantidad de complejidad en el modelo? Hacemos esto mediante la creación de un modelo en el que p = c, que es la probabilidad de capturar un animal que ya está marcado no es diferente a la captura de un animal sin identificación. Esto requiere de 2 líneas: una para especificar una fórmula para saber cómo se modela p, y la segunda que incorpora esa fórmula en una lista (en realidad una "lista de listas"). Utilizamos la 'share = TRUE "convención porque esto permitirá que p y c (técnicamente diferentes tipos de parámetros) para tomar en valores comunes. La fórmula ~ 1 especifica una intersección única para la matriz de diseño, por lo constantes p = c.
>pdotshared=list(formula=~1,share=TRUE)
> m0<-mark(edwards.eberhardt,model="Closed", model.parameters=list(p=pdotshared))
STOP NORMAL EXIT
....
> m0$results$real
estimate se lcl ucl fixed note
p g1 t1 0.081955 0.008816 0.0662522 0.1009771
N g1 96.258772 6.878606 86.6032790 114.7066900
> m0$results$AICc
[1] 379.6029
Tenga en cuenta que el valor estimado de N es un poco más alto (alrededor de 96), pero también que el valor de AIC es menor, 379,6 vs 381,5. Como hemos visto en sesiones anteriores, utilizamos AIC para guiar la selección del modelo o, preferiblemente, podemos simplemente promedio a través de los modelos alternativos, y tomar en cuenta la incertidumbre del modelo. Volveremos a esto después de que corremos algunos modelos más.
Vamos a ejecutar un modelo adicional, esta vez permitiendo p varíe entre ocasiones (tiempo), pero p = c contrario de otra manera
> ptimeshared=list(formula=~time,share=TRUE)
> mt<-mark(edwards.eberhardt,model="Closed",model.parameters=list(p=ptimeshared))
> mt$results$real
estimate se lcl ucl fixed note
p g1 t1 9.465890e-02 3.073600e-02 4.921960e-02 1.743547e-01
p g1 t2 8.414130e-02 2.906480e-02 4.202310e-02 1.613627e-01
p g1 t3 9.465900e-02 3.073600e-02 4.921970e-02 1.743549e-01
p g1 t4 1.472477e-01 3.775670e-02 8.740740e-02 2.373984e-01
p g1 t5 8.414130e-02 2.906480e-02 4.202310e-02 1.613626e-01
p g1 t6 5.258780e-02 2.318180e-02 2.181250e-02 1.213952e-01
p g1 t7 1.893189e-01 4.228050e-02 1.197932e-01 2.860811e-01
p g1 t8 1.156946e-01 3.377600e-02 6.410890e-02 1.999219e-01
p g1 t9 4.207010e-02 2.079490e-02 1.572280e-02 1.077360e-01
p g1 t10 3.155280e-02 1.806120e-02 1.012560e-02 9.401660e-02
p g1 t11 1.682832e-01 4.011370e-02 1.034393e-01 2.619019e-01
p g1 t12 5.258820e-02 2.318190e-02 2.181280e-02 1.213957e-01
p g1 t13 2.103520e-02 1.478950e-02 5.230800e-03 8.071680e-02
p g1 t14 7.362350e-02 2.726840e-02 3.502820e-02 1.482131e-01
p g1 t15 9.465910e-02 3.073610e-02 4.921970e-02 1.743551e-01
p g1 t16 2.937879e-09 5.694183e-06 -1.115766e-05 1.116354e-05
p g1 t17 4.207070e-02 2.079510e-02 1.572310e-02 1.077369e-01
p g1 t18 1.051769e-01 3.230150e-02 5.658920e-02 1.872039e-01
N g1 9.507803e+01 6.613440e+00 8.585798e+01 1.129215e+02
>
> mt$results$AICc
[1] 354.5968
Ahora tenemos una AIC inferior (el más bajo aún) y un valor de N entre los último 2 (95 vs 94 y 96). Tal vez nos estamos acercando.
Hasta ahora hemos llevado a cabo 3 de los modelos "clásicos" (o sus equivalentes) que se analizan por Otis et al. (1978):
Los resultados de AIC para estos modelos se resumen a continuación
model npar AICc DeltaAICc weight Deviance
3 p(~time)c()N(~1) 19 354.5968 0.00000 9.999949e-01 305.2648
1 p(~1)c()N(~1) 2 379.6029 25.00616 3.715168e-06 364.8259
2 p(~1)c(~1)N(~1) 3 381.5498 26.95302 1.403539e-06 364.7640
y podemos usar los pesos de AIC weights para calcular un estimido de promedio ponderado para N (codificación en R aquí).
Ahora, vamos a considerar el uso de AIC para la inferencia de modelos múltiples en más detalle. Como se ha dicho antes, hay básicemente dos tipos de aproximaciones sobre la inferencia de modelos múltiples. En el primero se usa AIC para ordenar los modelos y seleccionar el mejor modelo, o por lo menos concentrar en un grupo de modelos con peso de AIC grandes. Pero, de hecho no es necesario seleccionar ningún modelo, especialemente si el énfasis es en la estimación de parámetros, en este caso la abundancia N. Bajo esta aproximación no se usa AIC para seleccionar un o más modelos, en lugar se usa AIC ara computar un valor estimado de promedio ponderado, dónde los pesos son los pesos de AIC. Esta aproximación usa los estimadoros bajo todos los modelos, y computa un valor estimado de promedio ponderado y su varianza no condicional (es decir, no depende en supuesto de un modelo específico).
Se puede collecionar los resultados de todos los modelos y relativemente fácilmente computar los pesos de AIC y el valor estimado de promedio ponderado y la varianza de N. Pero estos cálculos son implemtados (más o menos) automaticamente en RMark y MARK. Los pasos en RMark son relativemente fácil-
Crear un objeto de resultados de los models por el comando results<-collect.models(). Esta función colleciona todos los modelos que se han sido ejecutados, por eso es importanted que solo los models de un típo de estructura de datos están incluido. Si hay duda, borrar la memoria y reejecutar los modelos.
Aplicar el comando model.average(results) al objeto de resulados y guardarlo en un otro objecto, por ejemplo
avg<-model.average(results)
Obtenir los valores estimados de los parámetros reales (los p, c, y N) y sus varianzas del objeto. En el caso de N es un poco más complejo porque el valor estimado en el objeto no es de N pero de los animals no marcados, es decir N - M dónde M es el número de animales no marcados. Pero este valor es siempre conocido en estudios de CMR para poblaciones cerradas.
Estos paso y todo lo anterior están illustrado en el código de R para el ejemplo.
Note que los modelos son implementados dentro una function run.edwards.eberhardt(). Esta aproximación no es necesaria, es sola para conveniencia. Por ejemplo, sería fácil modificar la función para hacerla más general en entrado, p. ej. introduciendo el argumeto data.input en la funcíon
run.edwards.eberhardt(data.input)
y sustituyendo data.input por edward.eberhardt dentro la función. En esta manera se puede correr los mismos modelos para diferentes conjuntos de datos, solo, p.ej.
run.edwards.eberhardt(edwards.eberhardt)
run.edwards.eberhardt(otros.datos)
Por supuesto, podria ser más simple solo usar los comandos individuales, sin una función, p.ej. los comandos
>data(edwards.eberhardt)
>pdotshared=list(formula=~1,share=TRUE)
>ptimeshared=list(formula=~time,share=TRUE)
>m0<-mark(edwards.eberhardt,model="Closed",model.parameters=list(p=pdotshared))
>mb<-mark(edwards.eberhardt,model="Closed")
>mt<-mark(edwards.eberhardt,model="Closed",model.parameters=list(p=ptimeshared))
ejecutarán los modelos M0,Mb, y Mt.
Los pasos en el MARK GUI para implementar promedios sobre los models son diferentes y son implemtados usando el menú. Proporciono una comparación entre MARK y RMark para estos ejemplos más adelante.
Véase también el Apéndice 3 del libro de MARK en línea
Siguiente: Diseño de estudios