Como se nota en la conferencia, la ocupación de múltiples temporada depende de los datos recogidos en el marco del Diseño robusto, por lo que s sitios se muestrean en cada una de las T ocasiones primarios y J secundarias para J * T total de muestras por sitio (en J en general puede variar entre primaria ocasiones, pero por ahora nos asumen toda J son los mismos). Los parámetros necesarios para modelar estos datos son ahora:
psi (1) - la probabilidad de ocupación inicial al comienzo del estudio
ext (t) - probabilidad de que un sitio ocupado convierte desocupada entre períodos principales ty t +1
col (t) - la probabilidad de que un sitio desocupado queda ocupada entre períodos principales ty t +1
p (j, j = 1, J * T) - probabilidad de detección de un lado ocupado para cada una de las muestras J * T
Dados los primeros 3 parámetros de ocupación después del tiempo t = 1 se modela como
psi (t +1) = psi (t) * [1-ext (t)] + [1-psi (t)] * Col (t)
Los parámetros, a su vez, pueden ser modelados en términos de sitio y ocasión específica covariables o factores de agrupación.
Ejemplo de Grand Skink - modelado de ocupación de estaciones múltiples
Ilustro el modelado de ocupación la multi-estación con un ejemplo de Gran Eslizones, un reptil en peligro de extinción endémica de la región de Otago Central de la Isla Sur de Nueva Zelanda. Los datos se proporcionan como uno de los ejemplos de conjuntos de datos con la instalación del programa de PRESENCE desde el sitio PWRC, pero en su lugar se utiliza la función colext () en el paquete unmarked. Los datos se encuentran en un archivo CSV que contiene los registros de presencia / ausencia (columnas BP) durante 5 años ("estaciones") y con 3 estudios por temporada, en cada uno de 352 afloramientos rocosos ("tores"). Una sola covariable específica de sitio se proporciona en la columna de R, que es una variable que indica si el sitio está en pastos (1) o pastizales nativos (0).
Primero leemos los datos, especificar covariables a nivel de sitio y los efectos del tiempo en la temporada (anual) y los niveles de la encuesta:
> #install unmarked if not installed
> #install.packages("unmarked")
> #load library
> library(unmarked)
> #set here the directory where your data files are located
> data_dir<-"C:/Users/mike/Dropbox/company/pacblackduck/workshop/occupancy"
> #data_dir<-"C:/Documents and Settings/conroy/My Documents/Dropbox/company/pacblackduck/workshop/occupancy"
> #data_dir<-"C:/mydir"
> setwd(data_dir)
> data<-read.csv("Grand_Skinks.csv")
> y<-data[,2:16]
> S <- nrow(data) # number of sites
> J <- 3 # number of secondary sampling occasions
> T <- 5 # number of primary periods
> #site level covariates
> site.cov<-data.frame(pasture=data[,18])
> #dummy variables to create seasonal effects (need 2 for 3 occasions)
> t1<-rep(c(1,0,0,0,0),S)
> t2<-rep(c(0,1,0,0,0),S)
> t3<-rep(c(0,0,1,0,0),S)
> t4<-rep(c(0,0,0,1,0),S)
> yearly.cov<-data.frame(t1,t2,t3,t4)
> #additive time effects within season. to model effects change between seasons create an interaction term
> s1<-rep(c(1,0,0,1,0,0,1,0,0,1,0,0,1,0,0),S)
> s2<-rep(c(0,1,0,0,1,0,0,1,0,0,1,0,0,1,0),S)
> second.cov<-data.frame(s1,s2)
Estos se combinan en una trama de datos especializada:
> skink<-unmarkedMultFrame(y=y,siteCovs=site.cov, numPrimary=5,yearlySiteCovs=yearly.cov,obsCovs=second.cov)
>
Ya estamos preparados correr unos modelos:
> fm1<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~1, data=skink)
Warning message:
In colext.fit(formula = formula, data = data, J = J, method = method, :
NAs introduced by coercion
> #back transformations
> backTransform(fm1,'det')
Backtransformed linear combination(s) of Detection estimate(s)
Estimate SE LinComb (Intercept)
0.689 0.0231 0.796 1
Transformation: logistic
> backTransform(fm1,"psi")
Backtransformed linear combination(s) of Initial estimate(s)
Estimate SE LinComb (Intercept)
0.395 0.0315 -0.428 1
Transformation: logistic
> backTransform(fm1,"col")
Backtransformed linear combination(s) of Colonization estimate(s)
Estimate SE LinComb (Intercept)
0.0683 0.0119 -2.61 1
Transformation: logistic
> backTransform(fm1,"ext")
Backtransformed linear combination(s) of Extinction estimate(s)
Estimate SE LinComb (Intercept)
0.1 0.0183 -2.2
Transformation: logistic
Los mensajes de advertencia se limitan a indicar que algunos valores que faltan, es decir, algunos sitios no se tomaron muestras en muestras J * t. En general, nosotros nos encargaremos de los diseños de desequilibrio (como si hemos números de las encuestas entre estaciones variando) que permiten a estos valores para faltar. En nuestro ejemplo, tenemos un estimado de probabilidad inicial de ocupación de alrededor de 0.4, la colonización de aproximadamente 0.06, y la extinción de 0.1.
Vamos a proceder a la construcción de varios modelos más que permiten combinaciones de tiempo y efectos covariables sitio.
> #additional models
> #constant initial psi , gamma , and epsilon
> #secondary occasion specific detection
> fm2<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~s1+s2, data=skink)
Warning message:
In colext.fit(formula = formula, data = data, J = J, method = method, :
NAs introduced by coercion
>
> #constant initial psi , gamma , and epsilon
> #primary occasion specific detection
> fm3<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~t1+t2+t3+t4, data=skink)
Warning message:
In colext.fit(formula = formula, data = data, J = J, method = method, :
NAs introduced by coercion
>
> #constant initial psi , gamma , and epsilon
> #primary*secondary occasion specific detection
> fm4<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~(t1+t2+t3+t4)*(s1+s2), data=skink)
Warning message:
In colext.fit(formula = formula, data = data, J = J, method = method, :
NAs introduced by coercion
>
>
> #initial psi as function of pasture
> fm5<-colext(psiformula = ~pasture, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~1, data=skink)
Warning message:
In colext.fit(formula = formula, data = data, J = J, method = method, :
NAs introduced by coercion
>
> #constant initial psi , gamma and epsilon time dependent
> #
> fm6<-colext(psiformula = ~1, gammaformula = ~t1+t2+t3+t4, epsilonformula = ~t1+t2+t3+t4,
+ pformula =~1, data=skink)
Warning message:
In colext.fit(formula = formula, data = data, J = J, method = method, :
NAs introduced by coercion
>
> fm7<-colext(psiformula = ~pasture, gammaformula = ~t1+t2+t3+t4, epsilonformula = ~t1+t2+t3+t4,
+ pformula =~t1+t2+t3+t4+pasture, data=skink,control=list(maxit=1e4))
Warning message:
In colext.fit(control = list(maxit = 10000), formula = formula, :
NAs introduced by coercion
>
> fm8<-colext(psiformula = ~pasture, gammaformula = ~1, epsilonformula = ~1,
+ pformula =~t1+t2+t3+t4+pasture, data=skink,control=list(maxit=1e4))
Warning message:
In colext.fit(control = list(maxit = 10000), formula = formula, :
NAs introduced by coercion
>
> fms <- fitList('psi0(.)col(.)ext(.)p(.)' = fm1,'psi0(.)col(.)ext(.)p(s)'=fm2, 'psi0(.)col(.)ext(.)p(t)'=fm3,'psi0(.)col(.)ext(.)p(t*s)'=fm4,'psi0(past)col(.)ext(.)p(.)'=fm5,'psi0(.)col(t)ext(t)p(.)'=fm6,'psi0(past)col(t)ext(t)p(t+past)'=fm7,'psi0(past)col(.)ext(.)p(t+past)'=fm8)
> modSel(fms)
nPars AIC delta AICwt cumltvWt
psi0(past)col(.)ext(.)p(t+past) 10 1750.05 0.00 6.7e-01 0.67
psi0(past)col(.)ext(.)p(.) 5 1751.54 1.49 3.2e-01 0.98
psi0(past)col(t)ext(t)p(t+past) 18 1757.52 7.47 1.6e-02 1.00
psi0(.)col(.)ext(.)p(s) 6 1769.05 19.01 5.0e-05 1.00
psi0(.)col(.)ext(.)p(t*s) 18 1769.84 19.79 3.4e-05 1.00
psi0(.)col(.)ext(.)p(t) 8 1771.98 21.93 1.2e-05 1.00
psi0(.)col(.)ext(.)p(.) 4 1775.02 24.97 2.5e-06 1.00
psi0(.)col(t)ext(t)p(.) 12 1781.43 31.38 1.0e-07 1.00
>
Una vez más, los mensajes de advertencia, simplemente indicaron que hubo valores perdidos en los datos y no son motivo de alarma. Tenga en cuenta que durante los últimos 2 modelos incluí una opción 'maxit' - esto era debido a que estos 2 modelos fallaron en converger en el número predeterminado de iteraciones permitidas.
El procedimiento de selección de modelos indica que los 3 mejores modelos todos contienen un efecto sobre la ocupación inicial de los pastos. Si queremos podemos explorar aún más esta relación, según el modelo promedio de este parámetro (o, mejor aún, las predicciones para los niveles especificados de pastos = 1 o 0). El procedimiento para hacerlo es idéntica a la que hicimos para la ocupación de una sola temporada.
Aquí están los datos y la escritura de R para este ejemplo.
Siguiente - Tamaños de muestreo