Práctica: Modelos de ocupación con el paquete unmarked
En este ejemplo se modela la ocupación por un especie de ardilla (Sciurus vulgaris) en Suiza.
Este ejemplo se basa en Pp. 590-600 del Capítulo 10 libro Kéry y Royle (2016). Además de traducirlo al español, hice algunos cambios menores en el código.
Para este ejemplo se emplea la función occ del paquete unmarked.
El proceso ecológico, en este caso la ocupación, se modela como:
\({z_i} \tilde {} Bernoulli(\psi)\)
Mientras que el proceso observacional, en este caso la detección, se modela como:
\({y_i{_j}} {|} {z_i} \tilde {} Bernoulli({z_i}p)\)
Instalar las siguientes librerias:
library(unmarked)## Loading required package: reshape## Warning: package 'reshape' was built under R version 3.2.5## Loading required package: lattice## Warning: package 'lattice' was built under R version 3.2.5## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 3.2.5library(raster)## Warning: package 'raster' was built under R version 3.2.5## Loading required package: sp## Warning: package 'sp' was built under R version 3.2.5##
## Attaching package: 'sp'## The following object is masked from 'package:unmarked':
##
## coordinates##
## Attaching package: 'raster'## The following objects are masked from 'package:unmarked':
##
## getData, projectionlibrary(rgdal)## Warning: package 'rgdal' was built under R version 3.2.5## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files:
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: (autodetected)
## WARNING: no proj_defs.dat in PROJ.4 shared files
## Linking to sp version: 1.2-3#library(AICcmodavg)PASO 1: Leer datos y preparar covariables
data <- read.table("SwissSquirrels.txt", header=TRUE)
head(data); tail(data)## X...spec.name coordx coordy ele route forest det071 det072 det073
## 1 Squirrel 922942 63276 450 6.2 3 0 0 0
## 2 Squirrel 928942 79276 450 5.5 21 0 0 0
## 3 Squirrel 928942 103276 1050 4.3 32 0 0 0
## 4 Squirrel 934942 95276 950 4.5 9 1 1 0
## 5 Squirrel 934942 111276 1150 5.4 35 0 0 0
## 6 Squirrel 946942 95276 550 3.6 2 0 0 0
## date071 date072 date073 dur071 dur072 dur073
## 1 20 46 54 230 300 280
## 2 19 50 56 180 180 180
## 3 22 55 61 175 195 175
## 4 22 49 70 300 255 240
## 5 15 43 60 180 215 210
## 6 21 40 49 200 200 195## X...spec.name coordx coordy ele route forest det071 det072 det073
## 260 Squirrel 1228942 79276 1850 5.5 82 0 1 0
## 261 Squirrel 1234942 71276 1850 5.4 81 1 0 0
## 262 Squirrel 1234942 127276 1850 4.7 82 1 1 1
## 263 Squirrel 1240942 135276 2450 4.4 0 0 0 NA
## 264 Squirrel 1246942 135276 1850 5.6 69 1 0 1
## 265 Squirrel 1252942 111276 1850 6.2 53 1 0 1
## date071 date072 date073 dur071 dur072 dur073
## 260 29 56 72 305 285 285
## 261 30 58 73 310 275 285
## 262 51 71 84 300 300 315
## 263 68 78 NA 282 304 NA
## 264 50 60 77 310 305 296
## 265 52 63 74 330 300 300View(data)
names(data)## [1] "X...spec.name" "coordx" "coordy" "ele"
## [5] "route" "forest" "det071" "det072"
## [9] "det073" "date071" "date072" "date073"
## [13] "dur071" "dur072" "dur073"str(data)## 'data.frame': 265 obs. of 15 variables:
## $ X...spec.name: Factor w/ 1 level "Squirrel": 1 1 1 1 1 1 1 1 1 1 ...
## $ coordx : int 922942 928942 928942 934942 934942 946942 946942 952942 958942 958942 ...
## $ coordy : int 63276 79276 103276 95276 111276 95276 111276 119276 111276 127276 ...
## $ ele : int 450 450 1050 950 1150 550 750 650 550 550 ...
## $ route : num 6.2 5.5 4.3 4.5 5.4 3.6 3.9 6.1 5.1 3.9 ...
## $ forest : int 3 21 32 9 35 2 6 60 5 13 ...
## $ det071 : int 0 0 0 1 0 0 0 0 0 0 ...
## $ det072 : int 0 0 0 1 0 0 0 0 0 0 ...
## $ det073 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ date071 : int 20 19 22 22 15 21 21 37 29 20 ...
## $ date072 : int 46 50 55 49 43 40 53 50 43 40 ...
## $ date073 : int 54 56 61 70 60 49 64 65 57 60 ...
## $ dur071 : int 230 180 175 300 180 200 150 210 205 180 ...
## $ dur072 : int 300 180 195 255 215 200 210 210 255 185 ...
## $ dur073 : int 280 180 175 240 210 195 135 190 180 195 ...Tabla 1. Datos de los muestreos de la ardilla en 265 sitios de muestreo, y covariables medidas. Aquí Solo se muestran los primeros 20 datos.
| X…spec.name | coordx | coordy | ele | route | forest | det071 | det072 | det073 | date071 | date072 | date073 | dur071 | dur072 | dur073 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Squirrel | 922942 | 63276 | 450 | 6.2 | 3 | 0 | 0 | 0 | 20 | 46 | 54 | 230 | 300 | 280 |
| Squirrel | 928942 | 79276 | 450 | 5.5 | 21 | 0 | 0 | 0 | 19 | 50 | 56 | 180 | 180 | 180 |
| Squirrel | 928942 | 103276 | 1050 | 4.3 | 32 | 0 | 0 | 0 | 22 | 55 | 61 | 175 | 195 | 175 |
| Squirrel | 934942 | 95276 | 950 | 4.5 | 9 | 1 | 1 | 0 | 22 | 49 | 70 | 300 | 255 | 240 |
| Squirrel | 934942 | 111276 | 1150 | 5.4 | 35 | 0 | 0 | 0 | 15 | 43 | 60 | 180 | 215 | 210 |
| Squirrel | 946942 | 95276 | 550 | 3.6 | 2 | 0 | 0 | 0 | 21 | 40 | 49 | 200 | 200 | 195 |
| Squirrel | 946942 | 111276 | 750 | 3.9 | 6 | 0 | 0 | 0 | 21 | 53 | 64 | 150 | 210 | 135 |
| Squirrel | 952942 | 119276 | 650 | 6.1 | 60 | 0 | 0 | 0 | 37 | 50 | 65 | 210 | 210 | 190 |
| Squirrel | 958942 | 111276 | 550 | 5.1 | 5 | 0 | 0 | 0 | 29 | 43 | 57 | 205 | 255 | 180 |
| Squirrel | 958942 | 127276 | 550 | 3.9 | 13 | 0 | 0 | 0 | 20 | 40 | 60 | 180 | 185 | 195 |
| Squirrel | 958942 | 143276 | 1150 | 3.8 | 50 | 0 | 0 | 0 | 14 | 40 | 52 | 230 | 220 | 195 |
| Squirrel | 964942 | 103276 | 750 | 7.7 | 57 | 0 | 1 | 1 | 21 | 57 | 70 | 240 | 225 | 240 |
| Squirrel | 964942 | 135276 | 1250 | 1.9 | 84 | 0 | 0 | 0 | 49 | 65 | 84 | 170 | 165 | 167 |
| Squirrel | 970942 | 95276 | 750 | 7.9 | 15 | 1 | 0 | 0 | 18 | 52 | 64 | 300 | 280 | 255 |
| Squirrel | 970942 | 127276 | 450 | 4.4 | 17 | 0 | 0 | 0 | 15 | 43 | 62 | 230 | 210 | 210 |
| Squirrel | 970942 | 159276 | 1050 | 6.0 | 58 | 1 | 0 | 0 | 16 | 43 | 56 | 135 | 160 | 145 |
| Squirrel | 976942 | 119276 | 750 | 6.2 | 26 | 0 | 1 | 0 | 26 | 49 | 69 | 230 | 230 | 220 |
| Squirrel | 976942 | 151276 | 1250 | 6.1 | 32 | 0 | 0 | 1 | 36 | 66 | 79 | 255 | 340 | 285 |
| Squirrel | 982942 | 63276 | 1250 | 4.5 | 66 | 1 | 1 | 0 | 24 | 53 | 72 | 255 | 255 | 240 |
| Squirrel | 982942 | 79276 | 350 | 6.2 | 45 | 1 | 0 | 0 | 22 | 49 | 61 | 180 | 195 | 210 |
# datos de presencia/ausencia
y <- as.matrix(data[,7:9]); head(y); tail(y)## det071 det072 det073
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 1 1 0
## [5,] 0 0 0
## [6,] 0 0 0## det071 det072 det073
## [260,] 0 1 0
## [261,] 1 0 0
## [262,] 1 1 1
## [263,] 0 0 NA
## [264,] 1 0 1
## [265,] 1 0 1# valores originales de las covariables sin estandarizar
elev.orig <- data[,"ele"]; head(elev.orig); tail(elev.orig)## [1] 450 450 1050 950 1150 550## [1] 1850 1850 1850 2450 1850 1850forest.orig <- data[,"forest"]; head(forest.orig); tail(forest.orig)## [1] 3 21 32 9 35 2## [1] 82 81 82 0 69 53time <- matrix(as.character(1:3), nrow=265, ncol=3, byrow=T); head(time); tail(time)## [,1] [,2] [,3]
## [1,] "1" "2" "3"
## [2,] "1" "2" "3"
## [3,] "1" "2" "3"
## [4,] "1" "2" "3"
## [5,] "1" "2" "3"
## [6,] "1" "2" "3"## [,1] [,2] [,3]
## [260,] "1" "2" "3"
## [261,] "1" "2" "3"
## [262,] "1" "2" "3"
## [263,] "1" "2" "3"
## [264,] "1" "2" "3"
## [265,] "1" "2" "3"date.orig <- as.matrix(data[,10:12]); head(date.orig); tail(date.orig)## date071 date072 date073
## [1,] 20 46 54
## [2,] 19 50 56
## [3,] 22 55 61
## [4,] 22 49 70
## [5,] 15 43 60
## [6,] 21 40 49## date071 date072 date073
## [260,] 29 56 72
## [261,] 30 58 73
## [262,] 51 71 84
## [263,] 68 78 NA
## [264,] 50 60 77
## [265,] 52 63 74dur.orig <- as.matrix(data[,13:15]);head(dur.orig); tail(dur.orig)## dur071 dur072 dur073
## [1,] 230 300 280
## [2,] 180 180 180
## [3,] 175 195 175
## [4,] 300 255 240
## [5,] 180 215 210
## [6,] 200 200 195## dur071 dur072 dur073
## [260,] 305 285 285
## [261,] 310 275 285
## [262,] 300 300 315
## [263,] 282 304 NA
## [264,] 310 305 296
## [265,] 330 300 300# histogramas
covs <- cbind(elev.orig, forest.orig, date.orig, dur.orig);head(covs); tail(covs)## elev.orig forest.orig date071 date072 date073 dur071 dur072 dur073
## [1,] 450 3 20 46 54 230 300 280
## [2,] 450 21 19 50 56 180 180 180
## [3,] 1050 32 22 55 61 175 195 175
## [4,] 950 9 22 49 70 300 255 240
## [5,] 1150 35 15 43 60 180 215 210
## [6,] 550 2 21 40 49 200 200 195## elev.orig forest.orig date071 date072 date073 dur071 dur072 dur073
## [260,] 1850 82 29 56 72 305 285 285
## [261,] 1850 81 30 58 73 310 275 285
## [262,] 1850 82 51 71 84 300 300 315
## [263,] 2450 0 68 78 NA 282 304 NA
## [264,] 1850 69 50 60 77 310 305 296
## [265,] 1850 53 52 63 74 330 300 300par(mfrow=c(2,2))
for(i in 1:8){
hist(covs[,i], breaks=50, col="grey", main=colnames(covs)[i])
}
pairs(cbind(elev.orig, forest.orig, date.orig, dur.orig))
———————————————–
PASO 2: estandarizacion de covariables
(means <- c(apply(cbind(elev.orig, forest.orig), 2, mean),
date.orig=mean(c(date.orig), na.rm=TRUE),
dur.orig=mean(c(dur.orig), na.rm=TRUE)))## elev.orig forest.orig date.orig dur.orig
## 1189.24528 34.49811 49.41232 231.87231(sds <- c(apply(cbind(elev.orig, forest.orig), 2, sd),
date.orig=sd(c(date.orig), na.rm=TRUE),
dur.orig=sd(c(dur.orig), na.rm=TRUE)))## elev.orig forest.orig date.orig dur.orig
## 639.81487 27.62343 21.58085 60.43774elev <- (elev.orig-means[1])/sds[1]; head(elev); tail(elev)## [1] -1.1554050 -1.1554050 -0.2176337 -0.3739289 -0.0613385 -0.9991098## [1] 1.032728 1.032728 1.032728 1.970499 1.032728 1.032728forest <- (forest.orig-means[2])/sds[2]; head(forest); tail(forest)## [1] -1.14026793 -0.48864723 -0.09043457 -0.92306103 0.01816888 -1.17646908## [1] 1.7196229 1.6834218 1.7196229 -1.2488714 1.2490080 0.6697896date <- (date.orig-means[3])/sds[3]; head(date); tail(date)## date071 date072 date073
## [1,] -1.362889 -0.15811774 0.21258121
## [2,] -1.409227 0.02723173 0.30525595
## [3,] -1.270215 0.25891858 0.53694280
## [4,] -1.270215 -0.01910564 0.95397912
## [5,] -1.594576 -0.29712985 0.49060543
## [6,] -1.316552 -0.43614196 -0.01910564## date071 date072 date073
## [260,] -0.94585303 0.3052560 1.046654
## [261,] -0.89951566 0.3979307 1.092991
## [262,] 0.07356910 1.0003165 1.602702
## [263,] 0.86130439 1.3246781 NA
## [264,] 0.02723173 0.4906054 1.278341
## [265,] 0.11990647 0.6296175 1.139329date[is.na(date)] <- 0; head(date); tail(date)## date071 date072 date073
## [1,] -1.362889 -0.15811774 0.21258121
## [2,] -1.409227 0.02723173 0.30525595
## [3,] -1.270215 0.25891858 0.53694280
## [4,] -1.270215 -0.01910564 0.95397912
## [5,] -1.594576 -0.29712985 0.49060543
## [6,] -1.316552 -0.43614196 -0.01910564## date071 date072 date073
## [260,] -0.94585303 0.3052560 1.046654
## [261,] -0.89951566 0.3979307 1.092991
## [262,] 0.07356910 1.0003165 1.602702
## [263,] 0.86130439 1.3246781 0.000000
## [264,] 0.02723173 0.4906054 1.278341
## [265,] 0.11990647 0.6296175 1.139329dur <- (dur.orig-means[4])/sds[4]; head(dur); tail(dur)## dur071 dur072 dur073
## [1,] -0.03097919 1.1272376 0.7963185
## [2,] -0.85827689 -0.8582769 -0.8582769
## [3,] -0.94100666 -0.6100876 -0.9410067
## [4,] 1.12723760 0.3826697 0.1344804
## [5,] -0.85827689 -0.2791685 -0.3618983
## [6,] -0.52735781 -0.5273578 -0.6100876## dur071 dur072 dur073
## [260,] 1.2099674 0.8790483 0.8790483
## [261,] 1.2926971 0.7135888 0.8790483
## [262,] 1.1272376 1.1272376 1.3754269
## [263,] 0.8294104 1.1934214 NA
## [264,] 1.2926971 1.2099674 1.0610538
## [265,] 1.6236162 1.1272376 1.1272376dur[is.na(dur)] <- 0; head(dur); tail(dur)## dur071 dur072 dur073
## [1,] -0.03097919 1.1272376 0.7963185
## [2,] -0.85827689 -0.8582769 -0.8582769
## [3,] -0.94100666 -0.6100876 -0.9410067
## [4,] 1.12723760 0.3826697 0.1344804
## [5,] -0.85827689 -0.2791685 -0.3618983
## [6,] -0.52735781 -0.5273578 -0.6100876## dur071 dur072 dur073
## [260,] 1.2099674 0.8790483 0.8790483
## [261,] 1.2926971 0.7135888 0.8790483
## [262,] 1.1272376 1.1272376 1.3754269
## [263,] 0.8294104 1.1934214 0.0000000
## [264,] 1.2926971 1.2099674 1.0610538
## [265,] 1.6236162 1.1272376 1.1272376———————————————–
PASO 3: procesar con unmarked
umf <- unmarkedFrameOccu(y=y,
siteCovs=data.frame(elev=elev, forest=forest),
obsCovs=list(time=time, date=date, dur=dur))
summary(umf)## unmarkedFrame Object
##
## 265 sites
## Maximum number of observations per site: 3
## Mean number of observations per site: 2.82
## Sites with at least one detection: 116
##
## Tabulation of y observations:
## 0 1 <NA>
## 553 194 48
##
## Site-level covariates:
## elev forest
## Min. :-1.46800 Min. :-1.24887
## 1st Qu.:-0.99911 1st Qu.:-0.95926
## Median :-0.06134 Median :-0.05423
## Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 1.03273 3rd Qu.: 0.77839
## Max. : 2.43938 Max. : 2.33504
##
## Observation-level covariates:
## time date dur
## 1:265 Min. :-1.6873 Min. :-2.4301
## 2:265 1st Qu.:-0.7605 1st Qu.:-0.6101
## 3:265 Median : 0.0000 Median : 0.0000
## Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.7223 3rd Qu.: 0.5481
## Max. : 3.1782 Max. : 5.5946# Modelos para la probabilidad de deteccion p(covariables) y psi(.)
summary(fm1 <- occu(~1 ~1, data=umf))##
## Call:
## occu(formula = ~1 ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.168 0.173 0.974 0.33
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -0.144 0.136 -1.06 0.291
##
## AIC: 805.9027
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 28
## Bootstrap iterations: 0summary(fm2 <- occu(~date ~1, data=umf))##
## Call:
## occu(formula = ~date ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.237 0.176 1.35 0.179
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.195 0.135 -1.45 0.146983
## date -0.386 0.112 -3.45 0.000559
##
## AIC: 795.9805
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 22
## Bootstrap iterations: 0summary(fm3 <- occu(~date+I(date^2) ~1, data=umf))##
## Call:
## occu(formula = ~date + I(date^2) ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.246 0.177 1.39 0.165
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1453 0.162 -0.898 0.369385
## date -0.3914 0.112 -3.497 0.000471
## I(date^2) -0.0582 0.104 -0.560 0.575372
##
## AIC: 797.6654
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 31
## Bootstrap iterations: 0summary(fm4 <- occu(~date+I(date^2)+I(date^3) ~1, data=umf))##
## Call:
## occu(formula = ~date + I(date^2) + I(date^3) ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.243 0.176 1.38 0.167
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.18737 0.164 -1.1391 0.255
## date -0.10885 0.232 -0.4700 0.638
## I(date^2) -0.00807 0.115 -0.0703 0.944
## I(date^3) -0.14500 0.107 -1.3543 0.176
##
## AIC: 797.6246
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 31
## Bootstrap iterations: 0summary(fm5 <- occu(~dur ~1, data=umf))##
## Call:
## occu(formula = ~dur ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.217 0.176 1.23 0.219
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.185 0.134 -1.38 0.16906
## dur 0.393 0.149 2.65 0.00813
##
## AIC: 801.2893
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 17
## Bootstrap iterations: 0summary(fm6 <- occu(~date+dur ~1, data=umf))##
## Call:
## occu(formula = ~date + dur ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.275 0.181 1.51 0.13
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.240 0.136 -1.76 0.0782
## date -0.379 0.113 -3.35 0.0008
## dur 0.368 0.153 2.40 0.0162
##
## AIC: 791.9682
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 32
## Bootstrap iterations: 0summary(fm7 <- occu(~date+I(date^2)+dur ~1, data=umf))##
## Call:
## occu(formula = ~date + I(date^2) + dur ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.283 0.183 1.55 0.122
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1971 0.163 -1.21 0.22728
## date -0.3865 0.114 -3.39 0.00069
## I(date^2) -0.0514 0.107 -0.48 0.63135
## dur 0.3629 0.153 2.37 0.01771
##
## AIC: 793.7368
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 30
## Bootstrap iterations: 0summary(fm8 <- occu(~date+I(date^2)+I(date^3)+dur ~1, data=umf))##
## Call:
## occu(formula = ~date + I(date^2) + I(date^3) + dur ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.281 0.182 1.54 0.122
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.230 0.165 -1.389 0.1649
## date -0.129 0.235 -0.548 0.5836
## I(date^2) -0.015 0.115 -0.130 0.8969
## I(date^3) -0.134 0.110 -1.222 0.2218
## dur 0.351 0.154 2.284 0.0224
##
## AIC: 794.115
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 39
## Bootstrap iterations: 0summary(fm9 <- occu(~dur+I(dur^2) ~1, data=umf))##
## Call:
## occu(formula = ~dur + I(dur^2) ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.292 0.187 1.56 0.118
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.0738 0.149 -0.496 0.6201
## dur 0.3737 0.134 2.788 0.0053
## I(dur^2) -0.2010 0.087 -2.309 0.0209
##
## AIC: 798.5154
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 31
## Bootstrap iterations: 0summary(fm10 <- occu(~date+dur+I(dur^2) ~1, data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.348 0.19 1.83 0.0677
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.131 0.1501 -0.876 0.381286
## date -0.376 0.1119 -3.361 0.000775
## dur 0.384 0.1350 2.843 0.004462
## I(dur^2) -0.190 0.0892 -2.126 0.033473
##
## AIC: 789.0856
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 32
## Bootstrap iterations: 0summary(fm11 <- occu(~date+I(date^2)+dur+I(dur^2) ~1, data=umf))##
## Call:
## occu(formula = ~date + I(date^2) + dur + I(dur^2) ~ 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.353 0.191 1.85 0.0647
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.0984 0.1737 -0.566 0.571250
## date -0.3812 0.1126 -3.384 0.000713
## I(date^2) -0.0407 0.1069 -0.381 0.703558
## dur 0.3839 0.1351 2.841 0.004501
## I(dur^2) -0.1876 0.0895 -2.096 0.036048
##
## AIC: 790.9401
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 34
## Bootstrap iterations: 0summary(fm12 <- occu(~date+I(date^2)+I(date^3)+dur+I(dur^2) ~1, data=umf))##
## Call:
## occu(formula = ~date + I(date^2) + I(date^3) + dur + I(dur^2) ~
## 1, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## 0.351 0.191 1.84 0.0659
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.12941 0.1759 -0.7358 0.4618
## date -0.14421 0.2353 -0.6128 0.5400
## I(date^2) -0.00911 0.1144 -0.0797 0.9365
## I(date^3) -0.12390 0.1102 -1.1243 0.2609
## dur 0.37799 0.1362 2.7762 0.0055
## I(dur^2) -0.18383 0.0896 -2.0522 0.0402
##
## AIC: 791.5865
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 44
## Bootstrap iterations: 0# para ordenar modelos de acuerdo al AIC mas bajo a mayor
fms <- fitList("p(.)psi(.)" = fm1,
"p(date)psi(.)" = fm2,
"p(date+date2)psi(.)" = fm3,
"p(date+date2+date3)psi(.)" = fm4,
"p(dur)psi(.)" = fm5,
"p(date+dur)psi(.)" = fm6,
"p(date+date2+dur)psi(.)" = fm7,
"p(date+date2+date3+dur)psi(.)" = fm8,
"p(dur+dur2)psi(.)" = fm9,
"p(date+dur+dur2)psi(.)" = fm10,
"p(date+date2+dur+dur2)psi(.)" = fm11,
"p(date+date2+date3+dur+dur2)psi(.)" = fm12)
(ms <- modSel(fms))## nPars AIC delta AICwt cumltvWt
## p(date+dur+dur2)psi(.) 5 789.09 0.00 0.4612 0.46
## p(date+date2+dur+dur2)psi(.) 6 790.94 1.85 0.1825 0.64
## p(date+date2+date3+dur+dur2)psi(.) 7 791.59 2.50 0.1321 0.78
## p(date+dur)psi(.) 4 791.97 2.88 0.1091 0.88
## p(date+date2+dur)psi(.) 5 793.74 4.65 0.0451 0.93
## p(date+date2+date3+dur)psi(.) 6 794.11 5.03 0.0373 0.97
## p(date)psi(.) 3 795.98 6.89 0.0147 0.98
## p(date+date2+date3)psi(.) 5 797.62 8.54 0.0065 0.99
## p(date+date2)psi(.) 4 797.67 8.58 0.0063 0.99
## p(dur+dur2)psi(.) 4 798.52 9.43 0.0041 1.00
## p(dur)psi(.) 3 801.29 12.20 0.0010 1.00
## p(.)psi(.) 2 805.90 16.82 0.0001 1.00cbind(fm1@AIC,fm2@AIC,fm3@AIC,fm4@AIC,fm5@AIC,fm6@AIC,fm7@AIC,fm8@AIC,
fm9@AIC,fm10@AIC,fm11@AIC,fm12@AIC) # mejor modelo 6 ## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 805.9027 795.9805 797.6654 797.6246 801.2893 791.9682 793.7368
## [,8] [,9] [,10] [,11] [,12]
## [1,] 794.115 798.5154 789.0856 790.9401 791.5865# modelos psi(covariable) y p(mejor modelo)
summary(fm13 <- occu(~date+dur+I(dur^2) ~elev, data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.287 0.190 1.51 0.131
## elev -0.343 0.217 -1.58 0.114
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.123 0.1522 -0.808 0.41936
## date -0.300 0.1218 -2.462 0.01381
## dur 0.435 0.1401 3.103 0.00191
## I(dur^2) -0.176 0.0916 -1.918 0.05512
##
## AIC: 788.7067
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 34
## Bootstrap iterations: 0summary(fm14 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2), data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2), data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 1.224 0.336 3.65 2.65e-04
## elev -0.203 0.192 -1.06 2.91e-01
## I(elev^2) -1.179 0.282 -4.18 2.91e-05
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.120 0.159 -0.754 0.4507
## date -0.235 0.116 -2.034 0.0419
## dur 0.375 0.146 2.562 0.0104
## I(dur^2) -0.103 0.126 -0.815 0.4150
##
## AIC: 767.4348
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 36
## Bootstrap iterations: 0summary(fm15 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+ I(elev^3), data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + I(elev^3),
## data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 1.2134 0.338 3.593 0.000327
## elev 0.0702 0.490 0.143 0.885935
## I(elev^2) -1.1462 0.293 -3.906 0.000094
## I(elev^3) -0.2049 0.340 -0.602 0.547287
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1196 0.160 -0.749 0.4538
## date -0.2369 0.116 -2.050 0.0403
## dur 0.3698 0.147 2.508 0.0121
## I(dur^2) -0.0996 0.129 -0.772 0.4402
##
## AIC: 769.0527
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 40
## Bootstrap iterations: 0cbind(fm13@AIC, fm14@AIC, fm15@AIC) # model 14 with elev2 best## [,1] [,2] [,3]
## [1,] 788.7067 767.4348 769.0527# Check effects of forest and interactions
summary(fm16 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+forest, data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + forest,
## data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.933 0.470 1.98 0.047264
## elev -0.542 0.265 -2.04 0.040959
## I(elev^2) -0.816 0.378 -2.16 0.030725
## forest 1.439 0.398 3.62 0.000298
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1987 0.159 -1.248 0.2119
## date -0.1867 0.113 -1.657 0.0975
## dur 0.3166 0.144 2.201 0.0277
## I(dur^2) -0.0533 0.131 -0.406 0.6850
##
## AIC: 738.0041
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 42
## Bootstrap iterations: 0summary(fm17 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+forest+I(forest^2), data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + forest +
## I(forest^2), data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.925 0.411 2.25 2.42e-02
## elev -0.495 0.249 -1.99 4.64e-02
## I(elev^2) -0.641 0.369 -1.74 8.26e-02
## forest 1.382 0.303 4.56 5.06e-06
## I(forest^2) -0.341 0.275 -1.24 2.15e-01
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1559 0.167 -0.936 0.3494
## date -0.2009 0.115 -1.750 0.0800
## dur 0.3165 0.148 2.136 0.0327
## I(dur^2) -0.0383 0.146 -0.263 0.7926
##
## AIC: 738.9185
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 47
## Bootstrap iterations: 0summary(fm18 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+forest+I(forest^2)+elev:forest, data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + forest +
## I(forest^2) + elev:forest, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.867 0.399 2.17 2.97e-02
## elev -0.477 0.246 -1.94 5.28e-02
## I(elev^2) -0.633 0.364 -1.74 8.20e-02
## forest 1.435 0.299 4.81 1.54e-06
## I(forest^2) -0.419 0.250 -1.67 9.43e-02
## elev:forest 0.309 0.266 1.16 2.44e-01
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1260 0.164 -0.769 0.4421
## date -0.2063 0.115 -1.798 0.0721
## dur 0.3061 0.149 2.058 0.0396
## I(dur^2) -0.0436 0.143 -0.306 0.7599
##
## AIC: 739.5644
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 46
## Bootstrap iterations: 0summary(fm19 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+forest+I(forest^2)+elev:forest+elev:I(forest^2), data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + forest +
## I(forest^2) + elev:forest + elev:I(forest^2), data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.853 0.399 2.140 3.24e-02
## elev -0.680 0.330 -2.060 3.94e-02
## I(elev^2) -0.653 0.372 -1.757 7.89e-02
## forest 1.383 0.295 4.696 2.65e-06
## I(forest^2) -0.325 0.285 -1.139 2.55e-01
## elev:forest 0.185 0.280 0.662 5.08e-01
## elev:I(forest^2) 0.257 0.272 0.944 3.45e-01
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.1396 0.167 -0.836 0.4030
## date -0.2017 0.115 -1.759 0.0786
## dur 0.3106 0.148 2.096 0.0361
## I(dur^2) -0.0366 0.145 -0.253 0.8006
##
## AIC: 740.6342
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 52
## Bootstrap iterations: 0summary(fm20 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+forest+I(forest^2)+elev:forest+elev:I(forest^2)+I(elev^2):forest, data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + forest +
## I(forest^2) + elev:forest + elev:I(forest^2) + I(elev^2):forest,
## data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.778 0.346 2.248 0.0246
## elev -0.771 0.362 -2.130 0.0332
## I(elev^2) -0.544 0.388 -1.402 0.1610
## forest 0.431 0.436 0.987 0.3235
## I(forest^2) -0.209 0.254 -0.823 0.4103
## elev:forest 0.197 0.371 0.532 0.5947
## elev:I(forest^2) 0.347 0.318 1.091 0.2753
## I(elev^2):forest 1.174 0.468 2.509 0.0121
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.0985 0.160 -0.616 0.5379
## date -0.2135 0.116 -1.845 0.0651
## dur 0.3242 0.147 2.205 0.0275
## I(dur^2) -0.0192 0.131 -0.147 0.8830
##
## AIC: 735.5701
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 48
## Bootstrap iterations: 0summary(fm21 <- occu(~date+dur+I(dur^2) ~elev+I(elev^2)+forest+I(forest^2)+elev:forest+elev:I(forest^2)+I(elev^2):forest+ I(elev^2):I(forest^2), data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) ~ elev + I(elev^2) + forest +
## I(forest^2) + elev:forest + elev:I(forest^2) + I(elev^2):forest +
## I(elev^2):I(forest^2), data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.953 0.390 2.441 0.0146
## elev -0.817 0.359 -2.278 0.0227
## I(elev^2) -0.907 0.488 -1.858 0.0632
## forest 0.642 0.465 1.380 0.1675
## I(forest^2) -0.488 0.338 -1.442 0.1494
## elev:forest 0.221 0.345 0.642 0.5209
## elev:I(forest^2) 0.435 0.332 1.311 0.1899
## I(elev^2):forest 0.991 0.483 2.053 0.0401
## I(elev^2):I(forest^2) 0.537 0.449 1.196 0.2319
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.0931 0.159 -0.587 0.5572
## date -0.2174 0.115 -1.883 0.0597
## dur 0.3192 0.147 2.165 0.0304
## I(dur^2) -0.0213 0.132 -0.161 0.8721
##
## AIC: 736.1344
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 53
## Bootstrap iterations: 0cbind(fm16@AIC, fm17@AIC, fm18@AIC, fm19@AIC, fm20@AIC) # fm20 is best## [,1] [,2] [,3] [,4] [,5]
## [1,] 738.0041 738.9185 739.5644 740.6342 735.5701# Check for some additional effects in detection
summary(fm22 <- occu(~date+dur+I(dur^2)+elev ~elev+I(elev^2)+forest+I(forest^2)+elev:forest+elev:I(forest^2)+I(elev^2):forest, data=umf))##
## Call:
## occu(formula = ~date + dur + I(dur^2) + elev ~ elev + I(elev^2) +
## forest + I(forest^2) + elev:forest + elev:I(forest^2) + I(elev^2):forest,
## data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.778 0.347 2.241 0.0251
## elev -0.768 0.393 -1.953 0.0508
## I(elev^2) -0.545 0.394 -1.383 0.1665
## forest 0.431 0.436 0.987 0.3235
## I(forest^2) -0.209 0.254 -0.822 0.4112
## elev:forest 0.199 0.378 0.526 0.5990
## elev:I(forest^2) 0.347 0.322 1.078 0.2812
## I(elev^2):forest 1.173 0.470 2.494 0.0126
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.09895 0.162 -0.6108 0.5413
## date -0.21301 0.120 -1.7695 0.0768
## dur 0.32438 0.148 2.1973 0.0280
## I(dur^2) -0.01935 0.131 -0.1478 0.8825
## elev -0.00308 0.191 -0.0161 0.9871
##
## AIC: 737.5698
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 56
## Bootstrap iterations: 0summary(fm23 <- occu(~dur+I(dur^2)+date*(elev+I(elev^2)) ~elev+I(elev^2)+forest+I(forest^2)+elev:forest+elev:I(forest^2)+I(elev^2):forest, data=umf))##
## Call:
## occu(formula = ~dur + I(dur^2) + date * (elev + I(elev^2)) ~
## elev + I(elev^2) + forest + I(forest^2) + elev:forest + elev:I(forest^2) +
## I(elev^2):forest, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.5429 0.347 1.5664 0.1173
## elev -1.1879 0.646 -1.8380 0.0661
## I(elev^2) 0.2739 0.766 0.3576 0.7206
## forest 0.4354 0.437 0.9962 0.3192
## I(forest^2) -0.2738 0.270 -1.0148 0.3102
## elev:forest -0.0184 0.489 -0.0377 0.9699
## elev:I(forest^2) 0.6090 0.453 1.3452 0.1786
## I(elev^2):forest 1.6080 0.707 2.2748 0.0229
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.1583 0.210 0.755 0.4505
## dur 0.3282 0.144 2.275 0.0229
## I(dur^2) 0.0311 0.122 0.254 0.7994
## date -0.3420 0.187 -1.827 0.0677
## elev -0.0356 0.193 -0.185 0.8536
## I(elev^2) -0.5687 0.286 -1.990 0.0466
## date:elev -0.0236 0.174 -0.136 0.8922
## date:I(elev^2) 0.1803 0.236 0.763 0.4452
##
## AIC: 739.2419
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 65
## Bootstrap iterations: 0summary(fm24 <- occu(~dur+I(dur^2)+date*(elev+I(elev^2))+forest ~elev+I(elev^2)+forest+I(forest^2)+elev:forest+elev:I(forest^2)+I(elev^2):forest, data=umf))##
## Call:
## occu(formula = ~dur + I(dur^2) + date * (elev + I(elev^2)) +
## forest ~ elev + I(elev^2) + forest + I(forest^2) + elev:forest +
## elev:I(forest^2) + I(elev^2):forest, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.58991 0.370 1.595 0.1108
## elev -1.14616 0.736 -1.558 0.1192
## I(elev^2) 0.46098 0.940 0.490 0.6239
## forest 0.30626 0.468 0.655 0.5125
## I(forest^2) -0.26618 0.279 -0.954 0.3400
## elev:forest -0.00815 0.511 -0.016 0.9873
## elev:I(forest^2) 0.64085 0.496 1.292 0.1965
## I(elev^2):forest 1.66510 0.785 2.121 0.0339
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.00355 0.231 -0.0154 0.9878
## dur 0.33321 0.143 2.3322 0.0197
## I(dur^2) 0.02299 0.119 0.1935 0.8465
## date -0.33494 0.187 -1.7921 0.0731
## elev -0.13665 0.204 -0.6704 0.5026
## I(elev^2) -0.61398 0.288 -2.1343 0.0328
## forest 0.26961 0.159 1.6948 0.0901
## date:elev -0.05251 0.173 -0.3035 0.7615
## date:I(elev^2) 0.16679 0.236 0.7060 0.4802
##
## AIC: 738.3106
## Number of sites: 265
## optim convergence code: 0
## optim iterations: 69
## Bootstrap iterations: 0cbind(fm22@AIC, fm23@AIC, fm24@AIC) # None better, hence, stay with model 20## [,1] [,2] [,3]
## [1,] 737.5698 739.2419 738.3106# -----------------------------
# OJO: modelo fm16
# para ver coeficientes
fm16[1]@estimates[1] #intercepto## (Intercept)
## 0.9332777fm16[1]@estimates[2] #elev## elev
## -0.5424226fm16[1]@estimates[3] #forest## I(elev^2)
## -0.816405# para linearizar psi y p
#backTransform(linearComb(fm16, coefficients = c(1,0,0), type= "state"))
#backTransform(linearComb(fm16, coefficients = c(1,0,0), type= "det"))———————————————–
PASO 5: para graficar
# se crean primero nuevas covariables estandarizadas
orig.elev <- seq(200, 2500,,100)
orig.forest <- seq(0, 100,,100)
orig.date <- seq(15, 110,,100)
orig.duration <- seq(100, 550,,100)
# estandarizadas
ep <- (orig.elev - means[1]) / sds[1]
fp <- (orig.forest - means[2]) / sds[2]
dp <- (orig.date - means[3]) / sds[3]
durp <- (orig.duration - means[4]) / sds[4]
# prediccion por covariable separada
newData <- data.frame(elev=ep, forest=0)
pred.occ.elev <- predict(fm20, type="state", newdata=newData, appendData=TRUE)
newData <- data.frame(elev=0, forest=fp)
pred.occ.forest <- predict(fm20, type="state", newdata=newData, appendData=TRUE)
newData <- data.frame(date=dp, dur=0)
pred.det.date <- predict(fm20, type="det", newdata=newData, appendData=TRUE)
newData <- data.frame(date=0, dur=durp)
pred.det.dur <- predict(fm20, type="det", newdata=newData, appendData=TRUE)
# graficos
par(mfrow=c(2,2), mar=c(5,5,5,5), cex.lab=1.2)
plot(pred.occ.elev[[1]] ~ orig.elev, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. occupancy prob.", xlab="Elevation (m)", frame=F)
matlines(orig.elev, pred.occ.elev[,3:4], lty=1, lwd=1, col="grey")
plot(pred.occ.forest[[1]] ~ orig.forest, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. occupancy prob.", xlab="Forest cover (%)", frame=F)
matlines(orig.forest, pred.occ.forest[,3:4], lty=1, lwd=1, col="grey")
plot(pred.det.date[[1]] ~ orig.date, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. detection prob.", xlab="Date (1= 1 April)", frame=F)
matlines(orig.date, pred.det.date[,3:4], lty=1, lwd=1, col="grey")
plot(pred.det.dur[[1]] ~ orig.duration, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. detection prob.", xlab="Survey duration (min)", frame=F)
matlines(orig.duration, pred.det.dur[,3:4], lty=1, lwd=1, col="grey")
# prediccion por covariables simultaneas
# Ocupacion ~ (forest, elevation) y deteccion ~ (survey duration, date)
pred.matrix1 <- pred.matrix2 <- array(NA, dim=c(100,100))
for(i in 1:100){
for(j in 1:100){
newData1 <- data.frame(elev=ep[i], forest=fp[j]) # para ocupacion
pred <- predict(fm20, type="state", newdata=newData1)
pred.matrix1[i, j] <- pred$Predicted
newData2 <- data.frame(dur=durp[i], date=dp[j]) # para deteccion
pred <- predict(fm20, type="det", newdata=newData2)
pred.matrix2[i, j] <- pred$Predicted
}
}
par(mfrow=c(1,1), mar=c(5,5,5,5), cex.lab=1.2)
mapPalette <- colorRampPalette(c("grey", "yellow", "orange", "red"))
image(x=orig.elev, y=orig.forest, z=pred.matrix1, col=mapPalette(100), axes=FALSE, xlab = "Elevation [m]", ylab="Forest cover [%]")
contour(x=orig.elev, y=orig.forest, z=pred.matrix1, add=TRUE, lwd=1.5, col="blue", labcex=1.3)
axis(1, at=seq(min(orig.elev), max(orig.elev), by=250))
axis(2, at=seq(0,100, by=10))
box()
title(main="Expected squirrel occurrence prob.", font.main=1)
points(data$ele, data$forest, pch="+", cex=1)
image(x=orig.duration, y=orig.date, z=pred.matrix2, col=mapPalette(100), axes=FALSE, xlab="Survey duration [min]", ylab="Date (1= April 1)")
contour(x=orig.duration, y=orig.date, z=pred.matrix2, add=TRUE, lwd=1.5, col="blue", labcex=1.3)
axis(1, at=seq(min(orig.duration), max(orig.duration), by=50))
axis(2, at=seq(0,100, by=10))
box()
title(main="Expected squirrel detection prob.", font.main=1)
matpoints(as.matrix(data[, 13:15]), as.matrix(data[, 10:12]), pch="+", cex=1)
——————————————————-
PASO 6: para mapear las ocupación
CH <- read.csv("Switzerland.csv", header=T)
View(CH)
# genera predicciones de la ocupación para cada cuadrante de 1 km2
newData <- data.frame(elev=(CH$elevation-means[1])/sds[1], forest=(CH$forest-means[2])/sds[2])
predCH <- predict(fm20, type="state", newdata=newData)## doing row 1000 of 42275
## doing row 2000 of 42275
## doing row 3000 of 42275
## doing row 4000 of 42275
## doing row 5000 of 42275
## doing row 6000 of 42275
## doing row 7000 of 42275
## doing row 8000 of 42275
## doing row 9000 of 42275
## doing row 10000 of 42275
## doing row 11000 of 42275
## doing row 12000 of 42275
## doing row 13000 of 42275
## doing row 14000 of 42275
## doing row 15000 of 42275
## doing row 16000 of 42275
## doing row 17000 of 42275
## doing row 18000 of 42275
## doing row 19000 of 42275
## doing row 20000 of 42275
## doing row 21000 of 42275
## doing row 22000 of 42275
## doing row 23000 of 42275
## doing row 24000 of 42275
## doing row 25000 of 42275
## doing row 26000 of 42275
## doing row 27000 of 42275
## doing row 28000 of 42275
## doing row 29000 of 42275
## doing row 30000 of 42275
## doing row 31000 of 42275
## doing row 32000 of 42275
## doing row 33000 of 42275
## doing row 34000 of 42275
## doing row 35000 of 42275
## doing row 36000 of 42275
## doing row 37000 of 42275
## doing row 38000 of 42275
## doing row 39000 of 42275
## doing row 40000 of 42275
## doing row 41000 of 42275
## doing row 42000 of 42275# para convertir a raster
PARAM <- data.frame(x=CH$x, y=CH$y, z=predCH$Predicted)
r1 <- rasterFromXYZ(PARAM)
# Máscara de cuadrantes con elevaciones mayores a 2250 m
elev <- rasterFromXYZ(cbind(CH$x, CH$y, CH$elevation))
elev[elev > 2250] <- NA
r1 <- mask(r1, elev)Tabla 2. Datos de los muestreos de los 42265 cuadrantes de 1 km2 de Suiza y covariables medidas. Aquí Solo se muestran los primeros 20 datos.
| X | x | y | elevation | forest | water |
|---|---|---|---|---|---|
| 1 | 910942 | 54276 | 340 | 6 | 0 |
| 2 | 910942 | 55276 | 340 | 11 | 1 |
| 3 | 911942 | 54276 | 380 | 4 | 0 |
| 4 | 911942 | 55276 | 390 | 72 | 3 |
| 5 | 911942 | 56276 | 357 | 9 | 7 |
| 6 | 911942 | 57276 | 357 | 5 | 5 |
| 7 | 911942 | 61276 | 500 | 0 | 0 |
| 8 | 911942 | 62276 | 472 | 43 | 0 |
| 9 | 911942 | 63276 | 462 | 6 | 0 |
| 10 | 911942 | 64276 | 472 | 0 | 0 |
| 11 | 912942 | 55276 | 401 | 30 | 0 |
| 12 | 912942 | 56276 | 362 | 12 | 0 |
| 13 | 912942 | 57276 | 367 | 4 | 0 |
| 14 | 912942 | 58276 | 379 | 5 | 7 |
| 15 | 912942 | 59276 | 382 | 12 | 9 |
| 16 | 912942 | 61276 | 483 | 0 | 0 |
| 17 | 912942 | 62276 | 464 | 30 | 0 |
| 18 | 912942 | 63276 | 456 | 14 | 0 |
| 19 | 912942 | 64276 | 450 | 4 | 0 |
| 20 | 913942 | 56276 | 390 | 2 | 1 |
———————————————————————
PASO 7:Estimación del área de ocurrencia de la ardilla durante el 2007
sum(predCH$Predicted) ## [1] 17354.57sum(predCH$Predicted[CH$elevation < 2250]) ## [1] 17350.34# estandariza covariables
pelev <- (CH$elevation-means[1])/sds[1]
pforest <- (CH$forest-means[2])/sds[2]
# función para predir la ocupación basado en el modelo 20
Eocc <- function(fm) {
betavec <- coef(fm)[1:8]
DM <- cbind(rep(1,length(pelev)), pelev, pelev^2, pforest, pforest^2, pelev*pforest, pelev*pforest^2, pelev^2*pforest)
pred <- plogis(DM%*%(betavec))
Eocc <- sum(pred)
Eocc
}
(Eocc.boot <- parboot(fm20, Eocc, nsim=1000, report=10))## t0 = 17354.57
## 17562,16936.3,18012.4,16407.6,19025.5,17282.5,19209.4,17087.2,19490.3,20408.7
## 17753.4,18859.3,18842.9,18885.5,16218.6,17507.9,17504.3,17515.9,16080.4,16630.5
## 18018.2,15165.1,14696.5,19005.4,19640.9,16498.1,18448.6,20185.5,17122.9,16905.4
## 17533,16736.7,16857.8,17989.1,17957.7,18203.2,17444.6,18778.5,18584.3,16408.8
## 16881.2,17402.1,16383,16363,18351.2,17555.9,16789.7,20181.6,19570.1,17466.3
## 18416.7,17930.4,18205.7,17579.3,19330.1,17695.1,18350.8,17920.6,16827.8,17853.9
## 18818.3,16178.4,16032.2,16787.4,15667.5,18904.5,17720,15192.2,18217.7,17525.7
## 14949.6,17413.3,18783,16662.6,20529.2,18389.5,19804.7,20874.8,16817.8,16552.9
## 15029.9,16681,16058.3,17200.9,17018.6,17294.8,16822.3,15460.4,17880.6,17397
## 18683.7,18028,17988.7,18292.5,19909.1,17307.6,15402.2,14808.9,18167,14887
## 16511.8,17597.3,17681.9,16979.8,18437.7,18299.7,15805.7,20244.8,13867,17098.7
## 16542.2,16307.3,19105.9,19542.8,17104.8,17423,16315.2,16536,16690.7,17251.9
## 17950.6,18515.1,17953.2,19144.7,16850,17302.1,16516.5,16843,20269.4,16356.1
## 18050.3,17761.6,18443.7,18803.2,18116.3,16959.8,17719.9,16171.2,17653.3,18986.4
## 16209,17304.6,17312.9,17406.9,17715.9,20803.6,15842.5,16540.9,16497.5,18890.6
## 18587.1,16567.8,17495,17164.9,16991.2,19256.1,17426.6,17310.9,15988.3,18678.9
## 17129.4,19309.3,17342.8,19503.2,17836.9,16195.6,15482.3,20690,16984.3,16361.5
## 16738.2,18944.2,18230.1,19821.2,16820.9,17138,19415.4,17955.8,17259.9,16727.6
## 16391.4,18123.4,17762,17764.7,17228,16580.2,15308.9,18689.2,19867.8,17399.9
## 19699.3,16541.2,16545.9,15325.8,20207.5,18929.6,16369.1,17905.7,20155.8,19805.9
## 15850.8,18142,17695.7,18066.5,18146.8,16096.7,16161.1,18541.9,16390.1,16740.3
## 17837.5,18168.6,18203.9,16829.4,17772.3,17197.4,16843.8,19502.3,17774.5,17077.6
## 17178.5,17749.2,16483.7,18056,19241.3,17794.6,18256,17666.1,21272.4,18229.8
## 16091.8,19868,16130.9,16391.1,19314.4,19173.6,16232.3,17146.8,17307,17362.7
## 16947.1,16624.4,15545.8,17154.4,17613.7,16471.9,17577.4,19847.3,15396.9,16298.5
## 17245.9,16661.8,15691.4,18119.6,17168.2,17371.4,19334.8,17901,16744.5,18135.5
## 17774.5,15971.3,16040.6,17358.8,16031,19793.2,18953.5,18018,18146.3,18215.6
## 16971,18591.7,16946.8,16854.7,17219.2,18019.2,15670,18776.4,18151.1,18680.5
## 17546.9,18798,16068.6,16261.8,17059.8,17568.6,17602.6,18629.9,20342.1,17984.5
## 16811.5,19167.5,14953.2,18927.4,17504.5,13874.2,18725.1,18482.3,18291.3,15847.2
## 18542.6,17016.3,17723.8,18334.9,18568.3,16248.2,16440,17985.1,17740.4,19209
## 17507.4,19415.9,18663.3,19129.6,18623,20622.6,17461.1,18467.5,17245.4,16592.6
## 18755.8,19847.5,16444.2,18425.7,17046.7,17178.8,17612.6,16821.9,20155.9,17953.2
## 14687.6,17391.4,16983.4,16877.2,15436.7,16466.6,18911.3,18974.2,15752.8,17981.3
## 18423.1,15482.5,17875.8,17363.2,16906.6,14612.7,18921.9,18185.7,15673.7,19127.4
## 18278.9,18126.1,17573.2,17394,18023.2,19816.7,14357.3,17337.9,18767,18761.1
## 18789.2,16955.4,19040.4,15905.9,16875.4,17339.6,17582.9,16160.5,16705.2,17456.7
## 17147.7,17726.5,16130,16684.6,18352.7,19560,16064.8,17920.3,18027.1,19128.3
## 15966.7,17112.5,17802.5,17936.3,19265.9,18053.5,16926.4,16111.7,17725.4,19462.9
## 17944.1,15657.2,16619.5,17904.7,16001.7,18130,17339.1,17506.1,18591,20139.1
## 18489.5,17199.2,17932.6,17888,15909.2,18306.4,17345.2,18820.8,17948.4,18751.6
## 16068.2,20344.4,16127.7,18528.3,16306.1,17147,15861.4,18931.5,16668,17954.5
## 15841.9,18905.6,16317.6,15571.1,18537.9,16204.3,17987,16073,16998.1,18229.3
## 17065.5,18139.9,17548.4,18077.9,18283,16120.7,17120.2,18197,18162.9,16391
## 19443.3,17760.5,19135.3,15134.2,17424.2,17866.8,17782.2,17851,17526.5,17362.1
## 16718,17776.1,17774.8,19231.7,16737.1,17232.3,17863.7,16345.8,21287.2,17045.3
## 16527.3,19294.4,19350.1,17013.9,16813.4,17243.4,15532.8,19008,18196.6,18540.4
## 15229.9,20638.1,17586.4,16077.6,15772.6,16500.7,18702,18326.5,16099.6,18090.1
## 17475.6,17781.4,17796.6,17082.9,17576.7,19016,15454.7,17508.4,19523.9,20378.4
## 17951.5,17130.3,17062.9,17378.3,15910,17136.2,17153,16102.4,17389.1,18571.2
## 18518.9,19128.4,19127.1,19055.4,17154.3,19243,14404.3,18621.5,18197.5,18113.5
## 16843.8,16532.6,16951.8,18230.9,17807.7,17487.4,18378,16352.8,18660.2,17869.2
## 17309.6,17227.1,18473,20390.6,15596.6,17277.6,17094.4,15399.4,19243.3,18566.7
## 16454.8,15660.8,14546.9,19532.6,18507.2,15145.5,16685.7,15263,18430.2,16527.5
## 16702.3,19270.6,18589.2,18602.8,18018.5,16967.3,18437,17674.7,19616.5,16861.9
## 17406.6,15640.7,17179.2,14985.8,18499.2,18758.8,17601.1,19397.4,15750,18327.3
## 18202.1,16038.2,18434.7,17165.5,16797.8,16521.7,17359.2,17312.3,17217,17492.6
## 19668.3,17400.6,15270.2,15723.5,15953.4,16613.9,17996.6,17297.6,16380.8,17837.1
## 15958.1,19053.4,18016.7,14370.9,18514.2,17551,17567.1,16136.2,15038.9,18354.5
## 15906.2,16472.1,18668.5,16627,17743,18963.1,16987.3,15640.3,18399.1,15291.6
## 16399.4,18002.6,17035.3,14968.9,19022.5,16283.6,18918.5,18487.5,19271.7,16964.2
## 16830.6,18196.1,17571.4,17852.6,17312.6,15540,14340.8,18531.2,18514.8,16846
## 20384.6,16738.6,16554.3,18492.2,19230.2,16359.9,16180.8,20661.3,16491.8,14132.7
## 19315.7,16681.1,16870.9,16407,14992.5,15830.4,17455.3,17274.4,17811,19034.1
## 17550.5,17653.9,16546.6,16384,17905.4,18341.8,16616.2,20752.3,19339.8,20239.3
## 18809.6,17047.1,18709.8,19419.3,17428.2,18383.9,15363.6,18289.2,16232.6,15341
## 17486.8,16896.4,19556.5,17795.6,17868,17907,18020.8,16526.6,15297.2,20970.4
## 18627.6,16347.3,16930.6,20972.5,18013.4,15868.7,19518,15675.2,17292,16901.3
## 20046.7,15540,19690.7,18133.9,18515.7,16996.4,18342.6,19532.5,19142.5,16844.2
## 16600.3,18749.4,17428.3,17847,18017.9,14771.5,16350,18161.2,16216.4,17232.1
## 17492.7,17725.1,16969.7,16005.9,18950.9,17022.6,17876.9,17473.1,18546.7,18614.1
## 18236.8,16630.6,18937.7,18015.5,18402.2,17545.8,17328.4,16234,17029.9,16347.5
## 16199.5,19939.9,19841.8,18529,18934.5,18086.7,18330.8,18291.6,16687.5,18494.6
## 16155.1,18262.9,17157,18384.5,17059.9,16497,15143.5,15782.5,15987,16397.8
## 19493.3,17928.6,17971.6,17172.1,15860.3,21049.5,18632.8,18198.5,16378.8,18175.6
## 16853.2,16506.8,18151.4,15061.7,16031.2,17751.4,17738.3,17067.5,17239.7,18447.5
## 17747.2,19664.4,17269.7,15165.5,16783.7,17469.1,16543.9,19673.3,17790.6,18203.4
## 17506.6,17435.4,17744.3,17367.4,17088.4,17447.6,17648.1,17445.3,18031.3,17015.1
## 18675.3,21934.4,20857.2,20329.4,18653.9,15797,18256,16307.6,19489.1,16600
## 17494.3,17780,15571.4,18644.9,16599.8,16947.8,18265.6,18912.6,17385.2,17664.4
## 15061.8,18116.2,16147.2,16037.3,16469.3,18850.4,14611.5,17664,16794.3,16730.5
## 15561.1,15028.3,14993.5,17100.6,19531.6,18197.3,16184.3,16658.2,16317.3,15336.8
## 17024.1,17748.1,16895.3,16452.8,16515.8,15883.2,15850.3,14951.3,18575.8,14393.8
## 17447.9,16566.5,17882.6,19010.9,15658.2,18899.9,16694.3,18869.4,17466.7,17063.9
## 16909.8,19505.2,18153.3,15809.6,19664.8,18376.4,16476.3,18409.9,17757.4,17435.3
## 13501.5,17904.6,17990.9,15926.6,15246.8,15777.3,17469.1,15807.7,18477.2,18008.4
## 17087.1,19043.2,17700.5,21488.5,17627.9,19149.2,18441.9,15205.7,17805.2,15778.5
## 16999.8,18505.3,18215,20287.4,16198.8,16275.6,17885.4,19736.4,18033.9,16120.9
## 17335.2,17740.7,16332.9,17182.4,17047.5,16112.9,18590,17267.6,17865.3,18855.9
## 18796,18386.8,16810.6,16690.6,14752,19187.4,16765.6,16659.7,17817.7,19041.5
## 17661,19592,16965.4,16883.9,16982.7,17337.4,15827.7,18602.9,17268.9,17073.7
## 16796.5,17787.6,19365.6,15985.1,15815.1,18168.7,16295.6,20176.9,18233.7,19043.9
## 19107.2,17763.4,17000.5,17547.3,17426.6,17460.8,16554.6,16948.1,17054.7,17129.5
## 18947.6,17465.5,18744.5,18142.5,16988.4,15548.9,16113.8,18220.2,14881.5,18357.9
## 16727.3,17869.8,18613.5,15803.3,18161.6,17425,17230.6,18922,19490.4,18854.5
## 16823.7,17178.9,17622.2,19072.4,18846,16336.2,17478.8,18640.8,18132.8,18494.9
## 17695.2,17286.6,16510.9,16854.8,19697.6,16824.8,17430.3,18158.6,17720.8,14988.7
## 17719.6,17073.6,16335.8,16487.1,17355.8,18422.1,16961.4,17587.5,16339.8,17123
## 15872.5,18081.7,16191.2,17353.6,19287.2,15939.8,16369.4,18268.9,17841,15380.2
## 20066.9,18821.7,18053.8,16540.7,16691.6,18852,16378.6,17403.6,19142.7,17101##
## Call: parboot(object = fm20, statistic = Eocc, nsim = 1000, report = 10)
##
## Parametric Bootstrap Statistics:
## t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)
## 1 17355 -178 1315 0.55
##
## t_B quantiles:
## 0% 2.5% 25% 50% 75% 97.5% 100%
## t*1 13501 14992 16619 17494 18404 20245 21934
##
## t0 = Original statistic compuated from data
## t_B = Vector of bootstrap samples—————————————————————-
RESULTADOS FINALES
# graficos (mismos de arriba)
par(mfrow=c(2,2), mar=c(5,5,5,5), cex.lab=1.2)
plot(pred.occ.elev[[1]] ~ orig.elev, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. occupancy prob.", xlab="Elevation (m)", frame=F)
matlines(orig.elev, pred.occ.elev[,3:4], lty=1, lwd=1, col="grey")
plot(pred.occ.forest[[1]] ~ orig.forest, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. occupancy prob.", xlab="Forest cover (%)", frame=F)
matlines(orig.forest, pred.occ.forest[,3:4], lty=1, lwd=1, col="grey")
plot(pred.det.date[[1]] ~ orig.date, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. detection prob.", xlab="Date (1= 1 April)", frame=F)
matlines(orig.date, pred.det.date[,3:4], lty=1, lwd=1, col="grey")
plot(pred.det.dur[[1]] ~ orig.duration, type="l", lwd=3, col="blue", ylim=c(0,1), las=1, ylab="Pred. detection prob.", xlab="Survey duration (min)", frame=F)
matlines(orig.duration, pred.det.dur[,3:4], lty=1, lwd=1, col="grey")
# Gráfico del bootstrap
plot(Eocc.boot)
# estimación área de ocurrencia en km2
(estimate.of.occurrence <- Eocc(fm20)) ## [1] 17354.57quantile(Eocc.boot@t.star, c(0.025, 0.975))## 2.5% 97.5%
## 14992.43 20245.42# proporción de área de ocurrencia
(c(point=17354.57, quantile(Eocc.boot@t.star, c(0.025,0.975)))/42275)## point 2.5% 97.5%
## 0.4105161 0.3546405 0.4788983# Mapa de distribución ocupación
par(mfrow=c(1,1), mar=c(3,3,3,3))
mapPalette <- colorRampPalette(c("grey","yellow","orange","red"))
plot(r1, col=mapPalette(100), axes=F, box=F, main="Distribucion de la ardilla en 2007")
points(data$coordx, data$coordy, pch="+", cex=0.8)
# Mapa distribución de incertidumbre
r2 <- rasterFromXYZ(data.frame(x=CH$x, y=CH$y, z=predCH$SE))
r2 <- mask(r2, elev)
plot(r2, col=mapPalette(100), axes=F, box=F, main="Incertidumbre de la distribucion")
points(data$coordx, data$coordy, pch="+", cex=0.8)
En conclusión, el análisis de ocupación de la ardilla en Suiza sugiere que esta especie ocupa 17,354 km2 lo cual corresponde al 41% de la superficie de este país. Las principales variables a la ocupación fueron la elevación (altitud) y la cobvertura del bosque.
Específicamente, a mayor altitud disminuye la ocupación; mientras que a coberturas mayores del 40% aumenta la ocupación. Por otro lado, los análisis sugieren que la probabilidad de detección, estimada en 0.xx, disminuye conforme pasan más días desde el inicio de muestreo, y aumenta conforme la duración de cada muestreo se incrementa.
Finalmente, es importante interpretar estos resultados en términos de la temporalidad. Los datos aquí analizados, aunque muy extensos, solo corresponden al muestreo de un solo año (2007). En este sentido, siempre es deseable y recomendable tener más años de muestreo para conocer variaciones en la ocupación, y definir si las mismas variables cada año afectan de similar manera la ocupación y la detección. Para estos casos, son convenientes los modelos de ocupación “multi-seasons”.