Configuración de un análisis de la ocupación en R

Aunque hay varias opciones para acceder a la modelización de ocupación desde dentro de R, el paquete unmarked proporciona una aproximación sencilla para el análisis de MLE. Además, los modelos de forma en que se formulan en unmarked son muy similares a la sintaxis de los métodos utilizados para la construcción de modelos lineales (lm) y modelos lineales generalizados (GLM) en R, por lo que ya debería estar familiarizado.

En primer lugar, asegúrese de que está ejecutando la versión más nueva de R. Si no es así, por favor, actualice a la última versión. En segundo lugar, porque unmarked no es parte de la instalación estándar de R, que tendrá que ir a uno de los espejos de CRAN e instalar la biblioteca antes de su puede construir modelos. Esto se hace desde el menú y la selección de paquetes \ Install Packages. Por otra parte, la línea de comandos R

> Install.packages ("unmarked")

se pueden introducir en la consola o desde un archivo de script.

Una vez instalado el paquete reside en el disco duro de su ordenador y no necesita volver a instalar (a menos que haya cambios), pero tendrá que ser cargado cada vez que R se ejecuta utilizando el mandato

> require (unmarked)

 o

> library (unmarked)

De nuevo, esto se puede introducir cada sesión desde una línea de comandos o archivo de script (más fácil), pero si se molesta usted para necesita cargar el paquete cada vez se inicializa el program, se puede modificar el perfil de R para cargarlo automaticamente.  Yo no hago esto de manera rutinaria porque por una cosa que no es una buena idea estar cargando un montón de bibliotecas en la memoria, a menos que realmente se va a estar utilizando en ese período de sesiones.

Una cierta cantidad de ayuda es disponible a través del comando integrado

>? unmarked

que proporciona la documentación en formato html. Un archivo pdf proporciona una cobertura un poco más detallada, y un artículo que acompaña una mayor elaboración.

Una vez que tienen los programas organizados, el primer paso para iniciar un análisis es tener los datos organizados de manera que sea legible por el programa. Un formato fácil de trabajar - es fácil porque es fácil de crear a partir de otros programas, como Excel, Access, etc, sino también porque es fácilmente digerible por desmarcado-es  un archivo de texto delimitado por comas  (*. xsv) que hemos estado trabajando con. Para un análisis de la ocupación, los datos de entrada caerán en 3 categorías:

• Historias de detección en la que las filas de los datos por lo general ser sitios y las columnas se representan ocasiones de detección u otras unidades de replicación,

• Covariables en sitios específicos, que son atributos medidos en cada uno de los sitios que potencialmente explican la variabilidad de sitio a sitio en la detección y / o la ocupación, y

• Específico de la muestra o covariables "nivel de observación", que son medidas de atributos en cada una de las ocasiones en las muestras (réplicas)

En los casos más simples (nuestros primeros ejemplos al menos) las covariables específicas del lugar no varían con el tiempo (por lo que se fijan en el tiempo para un sitio) y las covariables específicos de la muestra no varían entre los sitios, pero más generalmente covariantes sitio-específicas puede también varíar en ocasiones de las muestras y, a la inversa las covariables  específicas de la muestra pueden variar entre los sitios. Además, en nuestro ejemplo, vamos a leer las historias de detección y las covariables (por lo menos aquellos específicos del lugar) del mismo archivo de datos, pero es más fácil en muchos casos para almacenar estos diferentes tipos (y dimensiones) de datas en diferentes archivos.

Vamos a ilustrar la configuración de un análisis con los datos de ocupación del Blue Grosbeaks (caerulea Guiraca) en 41 antiguos campos plantados de pinos de hoja larga (Pinus palustris) en el sur de Georgia, EE.UU.. Las encuestas fueron 500 m transectos a través de cada campo y se completaron tres veces durante la temporada de cría en 2001. Los datos se encuentran en un archivo de texto delimitado por comas . La primera columna (que es opcional) es una etiqueta para cada sitio, simple numerada 1 - 41. Las próximas 3 columnas son las historias de detección para cada sitio en cada uno de 3 ocasiones durante la temporada de cría de 2001. Por ahora vamos a pasar por alto las covariables y otra información en la hoja de cálculo, y se centran en las historias de encuentros para cada sitio (columnas B, C y D a partir de la fila 3 y corriendo a la fila 43 (41 sitios en total).

Para leer los datos en el programa de R necesita saber dónde se encuentra en su computadora. Una manera sencilla de lograr esto es establecer un "directorio de trabajo" que albergará los datos y recibir archivos de entrada creados. Yo uso los comandos como éste - incluido en mi guión R:

>data_dir<-"C:/mydir"

>setwd(data_dir)

La primera línea crea un 'alias' llamada C :/ mydir, que puedo cambiar en cualquier momento. La segunda línea establece el directorio de trabajo con el nombre de alias.  Estoy primero voy a leer en sólo las historias de detección (columnas 2-4). Los comandos para hacer esto son:

>data<-read.csv("blgr.csv")

>#detection data rows are sites columns are detection replicates

>y<-data[,2:4]

>n<-nrow(data)

La línea penúltima crea un objeto que contiene sólo las historias de detección; es una matriz de 2 dimensiones con 41 filas y 3 columnas. La última línea es simplemente un contador de lo que puedo usar tarde para dimensionar otras partes de los datos (como covariables).

A continuación, vamos a leer en las covariables específicas del lugar en las columnas 5-9.

>#site level (individual) covariates

>blgr.site<-data[,5:9]

Esta es también una matriz de 2 dimensiones, pero con 41 filas y 4 columnas.

No hay covariables sitios específicos per se en este análisis, pero vamos a querer probar la variación específica de tiempo en la detección. Para ello, tendremos que crear una variable de tiempo que se comporta como un factor (en lugar de un valor numérico) en R. Podemos lograr esto mediante la creación de variables ficticias en una trama de datos, pero es más fácil definir un factor en R.

El código para hacer esto es simple:

> #create time factor and use as covariate

> #observation level (time specific) covariates

> time<-as.factor(rep(c(1,2,3),n))

> blgr.obs<-data.frame(time)

> 

Confirmamos que el tiempo es un factor:

> time

  [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1

 [38] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2

 [75] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

[112] 1 2 3 1 2 3 1 2 3 1 2 3

Levels: 1 2 3

Por último, ponemos las historias de detección y covariables en una sola trama de datos sin marcar llamado BLGR:

>#put everything together in unmarked data frame

>#note that covariate can come from separate files

>blgr <- unmarkedFrameOccu(y = y, siteCovs = blgr.site,obsCovs=blgr.obs)

>which can be summarized by

>#summary of unmarked data frame

>summary(blgr)

Este comando crea una tabla de estadísticos de resumen, siempre son útiles para asegurarse de que los datos se han leído correctamente (por ejemplo, debo esperar que todos los ceros o 1 de en el vector y, a menos que hubiera sitios no incluidos en la muestra algunas ocasiones (NA), pero no otros valores).

En realidad, esto no es una trama de datos estándar como hemos creado anteriormente, pero en vez de un tipo especial de R objeto llamado un objeto de tipo S4 que es directamente legible por el paquete no marcado como una trama de datos de ocupación. En este punto podríamos ahorrar este objeto y si queríamos tirar los datos originales (que no recomiendo eso!). Unos pocos puntos:

• Una vez más, podríamos haber leído las covariables desde diferentes archivos;

• Las covariables no se requieren si todo lo que estamos haciendo es ajustar el modelo de ocupación más simple (no hay efectos de sitio, la detección constante a lo largo de la muestra). En ese caso podríamos haber especificado 

  • BLGR <- unmarkedFrameOccu (y = y)

• El hecho de que tuvimos que crear covariables variables ficticias con el fin de modelar un efecto tiempo (que lo haremos más adelante) será molesto para la gente acostumbrada a los modelos pre-formados en MARK y presencia. Es simplemente el precio que se paga por la mayor flexibilidad - y similitud con GLM modelado - se obtiene en MARK.

• Por último, es posible que se haya preguntado acerca de las últimas 3 columnas de archivo de entrada csv, que no se han leído en MARK- unmarked. Estos son cargos de animales detectados en cada período de detección, y no se utilizan en el análisis básico de ocupación. Vamos a revisar estos si tenemos en cuenta los análisis que tienen que ver con la estimación de la abundancia a través de conteos repetidos.

Modelo básico: la detección y la ocupación homogénea

Los modelos de una sola estación que se están ejecutando en este punto se implementan en unmarked con la función occu(). Una vez que hemos establecido la trama de datos de ocupación en unmarked, es muy fácil de crear y ejecutar modelos. Para empezar, el modelo de detección y ocupación constante  (no hay tiempo, no hay efectos de covarianza) es creado por

>fm1<-occu(~1 ~1,blgr)

Esta sintaxis se puede utilizar una y otra vez, así que vamos a familiarizarse con él. En general, es

modelo <ocu-(~[formula para detección] ~ [formula para  ocupación],  trama de datos)

donde las fórmulas de detección y de ocupación son en forma de modelos lineales generalizados. En este caso, para ser explícitos estamos modelando la detección y la ocupación a través de modelos de intercepción (~ 1) utilizando una transformación logit, así por ejemplo, la probabilidad de detección

con una relación correspondiente, simple para Ψ probabilidad de ocupación

De hecho, ocu () proporciona estimadoes de estos 2 interceptiones 

lo que si (como suele ser el caso) queremos estimadoes de p y Ψ para realizar una logit inversa (expit) transformación. Unmarked ha construido en función backTranformation, que podemos usar aquí; más adelante nos encontraremos con que es más conveniente utilizar una función definida por el usuario que he construido. Volver al modelo:

>#CREATING MODELS

>#simple occupancy model

>fm1<-occu(~1 ~1,blgr)

when we type the command fm1

> fm1

Call:

occu(formula = ~1 ~ 1, data = blgr)

Occupancy:

 Estimate    SE    z P(>|z|)

     2.04 0.751 2.71 0.00666

Detection:

 Estimate    SE     z P(>|z|)

    0.206 0.242 0.851   0.395

AIC: 172.1898 

Así obtenemos estimados de los parámetros y un valor de AIC, pero se calculan los estimados en la escala logit. La función incorporada backTransform() proporcionará las estimaciones reales de escala de detección ('det') y ocupación ("estado")

> #back transformations

> backTransform(fm1,'det')

Backtransformed linear combination(s) of Detection estimate(s)

 Estimate     SE LinComb (Intercept)

    0.551 0.0599   0.206           1

Transformation: logistic 

> backTransform(fm1,"state")

Backtransformed linear combination(s) of Occupancy estimate(s)

 Estimate     SE LinComb (Intercept)

    0.885 0.0766    2.04           1

Transformation: logistic

Por último, unmarked proporciona un valor estimado emprico de Bayes e un intervalo de confianza del número / porcentaje de los sitios de estudio ocupados. Este método esencialmente utiliza la estimación de psi para generar un vector de variables binarias 0,1 z bajo la premisa de que z ~ Bernoulli (psi); la suma de la modalidad, más bajo, y cuantiles superiores de esta simulación proporciona entonces un CI empírica sobre Limitado / proporción de sitios ocupados. Para generar este durante los primeros estimaciones de los modelos que utilizamos

>

> 

> #empirical Bayes estimate of proportion of sites occupied

> s<-nrow(y)

> re <- ranef(fm1)

> EBUP <- bup(re, stat="mode")

> CI <- confint(re, level=0.95)

> print("95 % EB interval on number sites occupied")

[1] "95 % EB interval on number sites occupied"

> rbind(s_occup = c(Estimate = sum(EBUP), colSums(CI)) )

        Estimate 2.5% 97.5%

s_occup       33   33    41

> print("95 % EB interval on proportion of sites occupied")

[1] "95 % EB interval on proportion of sites occupied"

> rbind(PAO = c(Estimate = sum(EBUP), colSums(CI))/s )

    Estimate     2.5% 97.5%

PAO 0.804878 0.804878     1

Si se compara esta estadística para el valor estimado de back-tranformed de psi te darás cuenta de que es algo menor. Conceptualmente esto representa la proporción de la muestra finita (41 en este caso) de los sitios ocupados, mientras que psi estima la probabilidad de ocupación a partir de una lista infinita de sitios.

Adición de covariables para la detección y el modelado de ocupación

Teniendo en cuenta el marco de datos creado ya hemos establecido, es muy sencillo para añadir modelos contienen covariables para la predicción de la detección y de la ocupación. He creado 5 modelos más con diferentes combinaciones de detección y de ocupación covariables:

>#time specific detection, constant occupancy 

>fm2<-occu(~time ~1,blgr)

 >#constant detection, occupancy predicted by bqi 

>fm3<-occu(~1 ~bqi,blgr)

>#time specific detection, occupancy predicted by bqi

 >fm4<-occu(~time ~bqi,blgr)

>#constant detection, occupancy predicted by log(field size)

>fm5<-occu(~1 ~log(field.size),blgr)

>#constant detection, occupancy predicted by log(field size)+bqi

>fm6<-occu(~1 ~log(field.size)+bqi,blgr)

Podríamos construir más (si es que son una sensata priori) pero vamos a parar aquí.

En este punto vamos a querer resumir y comparar los modelos y tomar algunas decisiones. Un enfoque consiste en utilizar un procedimiento de selección de modelo integrado. La aplicación que aquí tenemos

>fmlist<-

fitList(p_dot_psi_dot=fm1,p_t_psi_t=fm2,p_dot_psi_bqi=fm3,p_t_psi_bqi=fm4,p_dot_psi_field=fm5,p_dot_psi

_field_bqi=fm6)

 > #model selection -- ranked list of models and AIC

> modSel(fmlist)

                    nPars    AIC delta AICwt cumltvWt

p_dot_psi_dot           2 172.19  0.00 0.381     0.38

p_dot_psi_bqi           3 173.12  0.93 0.239     0.62

p_dot_psi_field         3 173.92  1.73 0.160     0.78

p_dot_psi_field_bqi   4 175.11  2.92 0.088     0.87

p_t_psi_t                  4 175.30  3.11 0.081     0.95

p_t_psi_bqi               5 176.23  4.04 0.051     1.00

> 

Lo anterior indica que el modelo con el apoyo más alta es la p nula (.) Psi modelo (.). Usted puede ver estimados bajo distintos modelos escribiendo el nombre del modelo, por ejemplo

> fm3

Call:

occu(formula = ~1 ~ bqi, data = blgr)

Occupancy:

            Estimate   SE      z P(>|z|)

(Intercept)     2.69 1.41  1.909  0.0563

bqi            -1.39 1.55 -0.898  0.3690

Detection:

 Estimate    SE     z P(>|z|)

    0.206 0.242 0.851   0.395

AIC: 173.1211 

Sin embargo, estos valores estimados están en la escala logit. Puede obtener estimadoes a escala real a través de la función backTransform (), pero cuando covariables están presentes necesitará para especificar un valor para las covariables para su uso en la función lineal que predice p o psi. Lo que vamos a hacer en su lugar es tomar un enfoque de promedios modelo que tiene básicamente los siguientes pasos:

• Se especifica una cantidad que estamos tratando de predecir para cada modelo (por ejemplo, psi)

• Se especifica un valor o rango de valores de las covariables; menudo estos serán promedios sobre los datos, sino que puede ser cualquier valor.

• Generamos predicciones correspondientes a cada modelo y calculamos una media ponderada donde las ponderaciones son los pesos de la AIC.

• Calculamos el error estándar y el intervalo de confianza de la escala lineal

• Por último, transformamos los estimados y CI a la escala de probabilidad 

Vamos a seguir estos pasos con los datos de ejemplo y modelo establecido.

Comparación y inferencia de modelo múltiples

Unmarked tiene funciones incorporadas que hacen que la tarea de comparación de modelos y con un promedio de más fácil. Fundamentalmente, hay 3 pasos a seguir:

• Calcular AIC y pesas modelo para una lista específica de los modelos. Esto se hace por primera vez por especificando una lista de los modelos que ya han estado. Para nuestro ejemplo tenemos

>fmlist<-fitList(p_dot_psi_dot=fm1,p_t_psi_dot=fm2,p_dot_psi_bqi=fm3,p_t_psi_bqi=fm4,p_dot_psi_field=fm5,p_dot_psi_field_bqi=fm6)

•  A continuación,  pasar la lista a través de una función de selección de modelo

>modSel(fmlist)

• Finalmente especificamos una lista de valores de covarianza que queremos predecir las respuestas de, y el tipo de parámetro (ocupación, detección) que queremos predecir. Lo más sencillo es utilizar las covariables de todo el sitio se observan en los datos que produce la predicción

>#predict over observed covariates

> #predict over observed covariates

> new<-siteCovs(blgr)

> preddata.psi <- predict(fmlist, type = "state",new=new,appendData=T)

El preddata.psi trama de datos contiene ahora la ocupación prevista, se, y el intervalo de confianza para cada combinación de covariables sitio en los datos. Con frecuencia, sin embargo nos interesa sólo predicciones en ciertos niveles de las covariables observadas, no necesariamente en los datos. Por ejemplo podemos ver 4 combinación del efecto del tamaño del campo y BQI por:

> #create a new data frame to use for prediction for specified levels

> new<-data.frame(bqi=c(0,0,1,1),field.size=c(5,30,5,30))

> preddata.psi.new <- predict(fmlist, type = "state",new=new,appendData=T)

> preddata.psi.new

  Predicted         SE     lower     upper bqi field.size

1 0.9115003 0.08522441 0.5522803 0.9840954   0          5

2 0.8955429 0.09810999 0.5224137 0.9804994   0         30

3 0.8558248 0.12375452 0.5068063 0.9693401   1          5

4 0.8382170 0.12214459 0.5151838 0.9633650   1         30

Tenga en cuenta que entramos en el valor para el tamaño de campo en las unidades originales (ha) ya que la transformación de registro se lleva a cabo dentro del modelo.

 Todos los pasos anteriores se guardan en el archivo de script R adjunto, que se puede modificar para otras aplicaciones.

Siguiente: Tamaños del muestreo