Como acabamos de ver, tenemos unas aproximaciones de estimar el reclutamiento a partir de datos de CMR, pero éstas se ven obstaculizados tanto por el hecho de que los diseños de JS o Pradel y los modelos usando esta estructura de datos no tienen en cuenta la posible heterogeneidad de captura o efecto de comportamiento (captura previa) en recaptura. Estas influencias son conocidos causar graves sesgos en los valores estimados de abundancia o reclutamiento. El Diseño Robusto fue concebido originalmente para evitar este problema mediante un diseño en el que los datos se recogen de forma anidada:
Períodos primarios están espaciados a lo largo de intervalos durante los cuales se asume que la población sea abierto a adiciones (nacimientos, reclutamiento) y las pérdidas (muertes, emigración permanente).
Ocasiones secundarios suelen realizarse en períodos más cortos, en los que se asume la población sea cerrado.
La aplicación original del RD al problem de estimar reclutamiento luego utilizó modelos de CJS para estimar la supervivencia durante los períodos primarios, y modelos de CMR cerrados para estimar la abundancia N para cada período primario con las ocasiones secundarias como réplicas, y puso las 2 piezas de información en conjunto para estimar el reclutamiento. Las aproxaciones modernos realizan un análisis integrado de los dos tipos de datos. Vamos a ilustrar este enfoque mediante la aplicación de la extensión RD enfoque simetría temporal de Pradel.
En el ejemplo de los datos proviene de un estudio realizado por Jim Nichols, de sexo machos de Microtus pennsylvanicus atrapado en el centro de Patuxent, Maryland más de 6 meses, y los períodos de día durante 5 días cada mes (citado en los capítulos 17 a 19 de Williams et al.). Los datos están aquí, en formato de entrada en MARK.
Para analizar estos datos, hay que primero leerlos en R y dar nuevo formato para el programa RMark. Lo hago en el código adjunto, pero primero se analizan los datos secundarios dentro de cada período de primaria utilizando una serie de modelos cerrados. La razón para esto es tanto para hacer una verificación cruzada con el análisis robusto más complicado, pero también para proporcionar valores estimados preliminares de abundancia en caso de que estos son necesarios como valores iniciales.
Implementación en RMark
Pasando al análisis de diseño robusto, primero convertimos los datos de entrada a formato RMark
>#convert to RMark format
>mouse<-convert.inp(input.file)
A continuación, debemos especificar intervalos de tiempo para el Diseño Robusto. Aquí un "0" indica que todavía estamos en un período cerrado (así, a través de muestras secundarias dentro de una ocasión primaria), y un "1" indica que ha transcurrido un intervalo de tiempo (por ejemplo, un año) entre las muestras. Así:
Vamos a utilizar la parametrización de Pradel de Gamma (permanencia) de que es la forma original de este modelo. Como veremos más adelante, sin embargo, otras estadísticas quizá más útiles pueden ser el resultado de low valores estimados de Gamma.
>intervals<-c(rep(c(0,0,0,0,1),5),c(0,0,0,0))
>#set up Pradel model gamma parameterization
>mouseG.processed=process.data(mouse,model="RDPdGClosed",time.intervals=intervals)
mouseG.ddl=make.design.data(mouseG.processed)
A continuación, vamos a especificar algunas fuentes de variación en cada uno de los parámetros clave (p, c, Phi, y Gamma). Tenga en cuenta que por debajo de p y c (captura y recaptura probabilidades) "sesión" se refiere a los períodos entre ocasiones primarias, mientras que el "tiempo" se refiere a las ocasiones secundarias en las primarias.
Por otro lado el "tiempo" de Phi Gamma y siempre refieren las ocasiones primarias, ya que estos parámetros no tienen ningún significado dentro del período primario.
> #### seniority model model parameters
> #model parameters
> p=dot=list(formula=~1)
> cdot=list(formula=~1)
>
> p.t.s=dot=list(formula=~session+time)
> c.t.s=list(formula=~session+time)
> Phi.t=list(formula=~time)
> G.t=list(formula=~time)
> Phi.dot=list(formula=~1)
> G.dot=list(formula=~1)
>
A continuación, ponemos estos juntos en una serie de modelos; he corrido 5 aquí, pero muchos más son posibles. Tenga en cuenta que estos modelos son todos los modelos "cerrados" y no implican la heterogeneidad individual en la captura. Modelos correspondientes a los modelos de mezclas finitas que consideramos antes fácilmente se podrían construir. Por ejemplo especifando la estructura del modelo RDPdGFullHet lugar de RDPdGClosed en process.data () permitirá la heterogeneidad individual a través de la aproximación de mezcla.
Uno ejecutación de estos modelos podemos resumirlos por AIC y también obtener valores estimados del modelo promediada de los parámetros reales. Primero coleccionar los models en un objeto
> mods<-collect.models()
y después correr la función model.average para obtener estimados de promedio ponderado. Por ejemplo se proporcionan los valores de N por
mod_avg_N<-model.average(mods,"N")
Desde el objeto mods podemos extraer valores estimados de nuestros diversos parámetros de interés. Tenga en cuenta que para entender la indexación de los parámetros reales, es útil volver al objeto ddl..
Los valores estimados de promedio ponderado de abundancia, Phi, y Gamma se proporcionan por
> mod_avg_phi<model.average(mods,¨Phi¨)
y
> mod_avg_gamma<model.average(mods,¨Gamma¨) (notar las mayúsculas).
Tenga en cuenta que los valores estimados "reales" para N son realmente valores estimados de U: el número de animales no marcados en la población. Se obtiene N por añadiendo M, el número de animales marcados, a N = M + U. Por último, ya que éstas están cerradas (dentro del período primario) modelos, sabemos M simplemente del número de animales únicos que fueron capturados y marcados, que es sólo el número de filas de la matriz de la historia de la captura.
> mod_avg_N<-model.average(mods,"N")
>
> mod_avg_N<-subset(mod_avg_N,select=c(estimate,se))
> #THESE ARE UNMARKED MUST ADD MARKED ANIMALS EACH PERIOD
> ################
> mod_avg_N$marked<-marked
>
> #last 2 columns are estimates of N and SE
> mod_avg_N$Nhat<-marked+mod_avg_N$estimate
> mod_avg_N$Nhat.se<-mod_avg_N$se
>
Por último, los valores estimados del modelo promediado ponderado de Phi Gamma y también son capaces de producir valores estimados de la tasa de crecimiento de la población (Lambda) y el reclutamiento per cápita (f) por las relaciones fundamentales entre Gamma Phi, Lambda, y f resumen a continuación y elaborados en más detalle en Williams et al. (2002). Estas relaciones nos dan un instrumento simple de obtener valores estimados puntuales de Lambda y f dadas de Phi y Gamma; valores estimados de error estándar requieren métodos de aproximación tales como el método Delta, ya que por ejemplo el cálculo de lambda requiere división, una función no lineal.
>
> #model averaged gamma
> mod_avg_gamma<-estims[6:10,]
>
> #function to do delta approximation of variance of x/y
> delta_div<-function(est1,est2,se1,se2){
+ std.err<-array(dim=length(est1))
+ for(i in 1:length(est1))
+ {
+ fd<-c(1/est1[i],-est1[i]/est2[i]^2)
+ d<-matrix(data=fd,ncol=2)
+ v<-diag(c(se1[i]^2,se2[i]^2))
+ std.err[i]<-(sqrt(d%*%v%*%t(d)))
+ }
+ return(std.err)
+ }
>
> ests<-data.frame(phi=mod_avg_phi$estimate,se.phi=mod_avg_phi$se,gamma=mod_avg_gamma$estimate,se.gamma=mod_avg_gamma$se)
> ests$lambda<-ests$phi/ests$gamma
> ests$se.lambda<-delta_div(est1=ests$phi,est2=ests$gamma,se1=ests$se.phi,se2=ests$se.gamma)
> ests$f<-ests$lambda-ests$phi
> ests$se.f=sqrt(ests$se.lambda^2+ests$se.phi^2)
> print("Model averaged estimates of phi,gamma, f and lambda")
[1] "Model averaged estimates of phi,gamma, f and lambda"
> ests
phi se.phi gamma se.gamma lambda se.lambda f se.f
1 0.8626606 0.05151162 0.6486825 0.03011459 1.3298657 0.08589042 0.4672052 0.1001529
2 0.5547355 0.06245511 0.6518652 0.03274127 0.8509973 0.12042608 0.2962618 0.1356580
3 0.7288340 0.06799952 0.6496342 0.03042859 1.1219144 0.10708040 0.3930804 0.1268469
4 0.5905528 0.06232372 0.6502959 0.03119877 0.9081293 0.11417429 0.3175765 0.1300770
5 0.9497819 0.08284529 0.6466682 0.03231988 1.4687313 0.11400316 0.5189494 0.1409257
Aquí está el código R para realizar los análisis anteriores. La secuencia de comandos produce un archivo de entrada para su uso en la interfaz gráfica de usuario MARK, a continuación. Ofrezco esta información general, pero no es realmente necesario aquí, desde que empezamos con el mouse_robust.inp en el formato de MARK.
Implementación en MARK
Modelos correspondientes a los anteriores se pueden crear utilizando la interfaz gráfica de usuario MARK y el archivo mouse_robust.inp. El usuario debe seleccionar los datos y el tipo de modelo de la lista Data Type. Seleccione "Pradel Models Including Robust Designs" y (desde el panel siguiente) "Robust Design Pradel Survival and Seniority". Por último, seleccione "Full Likelihood p and c", lo que permitirá la modelización no mezcla de captura y recaptura similar a los modelos que hemos creado en R. A continuación, deberá especificar el número de los grupos (1) y las ocasiones de captura (30), y los intervalos de supervivencia sólidas apropiadas. Utilice la herramienta de Easy Robust Design (lado derecho de la pantalla) e introduzca 6 ocasiones primarios. El pop-up le dará la selección por defecto de 5 ocasiones secundarias cada primarias, que pueden ser retenidos para este ejemplo. A continuación, proceder al menú Run para seleccionar o construir modelos.
Estos pasos deberían estar familiarizados y no se repetirán aquí; yo sin embargo les he pedido para la construcción de algunos modelos en MARK como parte de los ejercicios de revisión.
Siguiente: Modelos de emigración temporales