Distance_Sampling
Oscar Ramírez-Alán
19 de abril de 2016
Se presenta una guía de cómo elaborar análisis de estimación de la abundancia, utilizando métodos de distancia (línea de transecto y puntos de conteo) y análisis de la covarianza para estimar esa variación que afecta la abundancia y densidad en los grupos de interés.
#Cargamos la base de datos que contiene la matrix principal de datos
dists.sub<-read.csv("Dista_line.csv", header=T, sep=";")
head(dists.sub,15)## date trans distance Cantidad
## 1 26-oct-12 L15 7.86 1
## 2 27-oct-12 L15 10.24 4
## 3 28-oct-12 L15 12.44 2
## 4 29-oct-12 L15 3.76 3
## 5 30-oct-12 L15 4.78 7
## 6 31-oct-12 L15 8.45 3
## 7 01-nov-12 L15 13.41 2
## 8 02-nov-12 L15 5.80 2
## 9 03-nov-12 L15 7.45 5
## 10 04-nov-12 L15 11.45 3
## 11 05-nov-12 L15 0.86 4
## 12 06-nov-12 L15 9.23 6
## 13 07-nov-12 L15 12.51 4
## 14 08-nov-12 L15 6.08 8
## 15 09-nov-12 L12 9.15 1dim(dists.sub) #nos permite ver la dimensión o confirmación de los datos## [1] 106 4dists.order <- dists.sub[order(dists.sub$dist), ] #Permite ordenas las distancias de menor a mayor (resulta util para establecer las truncaciones a utilizar)
dists.order # Observamos que la mayor distancia detectada en los datos es de 35.82 m## date trans distance Cantidad
## 78 11-ene-13 L17 0.02 1
## 79 12-ene-13 L17 0.60 1
## 58 22-dic-12 L14 0.62 3
## 55 19-dic-12 L14 0.75 1
## 53 17-dic-12 L14 0.81 1
## 11 05-nov-12 L15 0.86 4
## 102 04-feb-13 L124 0.97 2
## 46 10-dic-12 L14 0.99 3
## 61 25-dic-12 L95 0.99 2
## 62 26-dic-12 L95 1.00 2
## 63 27-dic-12 L95 1.22 4
## 85 18-ene-13 L17 1.26 2
## 54 18-dic-12 L14 1.49 1
## 96 29-ene-13 L17 1.71 3
## 86 19-ene-13 L17 1.72 4
## 28 22-nov-12 L14 1.85 3
## 67 31-dic-12 L95 1.91 2
## 80 13-ene-13 L17 2.00 3
## 47 11-dic-12 L14 2.28 2
## 30 24-nov-12 L14 2.61 2
## 41 05-dic-12 L11 2.91 3
## 43 07-dic-12 L14 2.92 2
## 97 30-ene-13 L17 3.08 4
## 42 06-dic-12 L14 3.25 1
## 68 01-ene-13 L95 3.26 1
## 37 01-dic-12 L13 3.31 3
## 98 31-ene-13 L17 3.31 3
## 25 19-nov-12 L16 3.58 5
## 44 08-dic-12 L14 3.69 1
## 69 02-ene-13 L95 3.73 4
## 4 29-oct-12 L15 3.76 3
## 40 04-dic-12 L11 3.77 4
## 50 14-dic-12 L14 3.82 3
## 18 12-nov-12 L16 3.85 4
## 33 27-nov-12 L14 4.01 2
## 24 18-nov-12 L16 4.21 2
## 60 24-dic-12 L14 4.37 4
## 64 28-dic-12 L95 4.61 1
## 20 14-nov-12 L16 4.70 2
## 72 05-ene-13 L95 4.75 4
## 5 30-oct-12 L15 4.78 7
## 105 07-feb-13 L124 4.86 2
## 36 30-nov-12 L13 5.10 4
## 23 17-nov-12 L16 5.11 4
## 77 10-ene-13 L95 5.79 4
## 8 02-nov-12 L15 5.80 2
## 70 03-ene-13 L95 5.83 2
## 71 04-ene-13 L95 5.90 3
## 38 02-dic-12 L11 5.97 2
## 14 08-nov-12 L15 6.08 8
## 31 25-nov-12 L14 6.15 4
## 92 25-ene-13 L17 6.26 1
## 16 10-nov-12 L12 6.38 2
## 103 05-feb-13 L124 6.62 3
## 35 29-nov-12 L13 6.90 3
## 81 14-ene-13 L17 6.94 3
## 82 15-ene-13 L17 7.21 1
## 9 03-nov-12 L15 7.45 5
## 59 23-dic-12 L14 7.61 3
## 74 07-ene-13 L95 7.64 3
## 83 16-ene-13 L17 7.66 1
## 1 26-oct-12 L15 7.86 1
## 87 20-ene-13 L17 8.39 3
## 6 31-oct-12 L15 8.45 3
## 15 09-nov-12 L12 9.15 1
## 12 06-nov-12 L15 9.23 6
## 65 29-dic-12 L95 9.25 3
## 34 28-nov-12 L14 9.73 5
## 32 26-nov-12 L14 9.74 3
## 93 26-ene-13 L17 10.02 3
## 57 21-dic-12 L14 10.03 2
## 84 17-ene-13 L17 10.19 4
## 56 20-dic-12 L14 10.21 2
## 2 27-oct-12 L15 10.24 4
## 75 08-ene-13 L95 10.59 2
## 52 16-dic-12 L14 11.08 3
## 26 20-nov-12 L14 11.16 4
## 10 04-nov-12 L15 11.45 3
## 27 21-nov-12 L14 12.23 2
## 73 06-ene-13 L95 12.37 4
## 94 27-ene-13 L17 12.39 3
## 3 28-oct-12 L15 12.44 2
## 104 06-feb-13 L124 12.45 2
## 13 07-nov-12 L15 12.51 4
## 19 13-nov-12 L16 12.55 4
## 90 23-ene-13 L17 12.78 1
## 45 09-dic-12 L14 13.21 3
## 91 24-ene-13 L17 13.24 1
## 88 21-ene-13 L17 13.35 3
## 48 12-dic-12 L14 13.38 3
## 7 01-nov-12 L15 13.41 2
## 22 16-nov-12 L16 14.55 1
## 106 08-feb-13 L124 15.41 4
## 66 30-dic-12 L95 15.81 4
## 49 13-dic-12 L14 16.23 2
## 100 02-feb-13 L17 16.64 4
## 76 09-ene-13 L95 17.83 8
## 21 15-nov-12 L16 17.86 3
## 39 03-dic-12 L11 18.37 2
## 51 15-dic-12 L14 19.28 4
## 99 01-feb-13 L17 19.42 3
## 101 03-feb-13 L13 19.42 1
## 89 22-ene-13 L17 19.44 2
## 95 28-ene-13 L17 19.46 1
## 17 11-nov-12 L12 21.21 3
## 29 23-nov-12 L14 35.82 4dim(dists.order)## [1] 106 4#Haremos una truncación a priori para utilizar las distancias menores a 20 metros
dists.trunc <- dists.order[which(dists.order$dist < 25), ]
dim(dists.trunc) ## [1] 105 4library(unmarked)## Warning: package 'unmarked' was built under R version 3.2.4## Loading required package: reshape
## Loading required package: lattice
## Loading required package: Rcpp## Warning: multiple methods tables found for 'overlay'library(Distance)## Loading required package: mrds
## This is mrds 2.1.15
## Built: R 3.2.2; ; 2016-04-01 21:22:19 UTC; windowslibrary(mrds)
yDat <- formatDistData(dists.trunc, distCol="distance", transectNameCol="trans", dist.breaks=c(0,5,10,15,20)) #Crea categoria de distancias (m). Esto cambia de acuerdo a su estudio.
yDat## [0,5] (5,10] (10,15] (15,20]
## L11 2 1 0 1
## L12 0 2 0 0
## L124 2 1 1 1
## L13 1 2 0 1
## L14 14 4 7 2
## L15 3 6 5 0
## L16 4 1 2 1
## L17 8 5 6 4
## L95 8 5 2 2nrow(yDat)## [1] 9hist(yDat)
Covariables
covs<-read.csv("Dista_line_cov.csv", header=T, sep=";")
covs## trans long type length hab area
## 1 L11 50 T 25 OW 2
## 2 L12 50 T 25 OW 1
## 3 L124 50 T 25 OW 2
## 4 L13 50 T 25 OW 1
## 5 L14 50 T 25 OW 2
## 6 L15 50 S 25 OG 1
## 7 L16 50 S 25 OG 2
## 8 L17 50 S 25 OG 1
## 9 L95 50 S 25 OG 1nrow(covs)## [1] 9#Unir las matrices de datos con las covariables
#Aqui es importante que no debe faltar ningun comando, si no la unión y análisis no podrá ejecutarse
umf <- unmarkedFrameDS(y=as.matrix(yDat), siteCovs=covs, survey="line",dist.breaks=c(0,5,10,15,20), tlength=covs$length, unitsIn="m")
umf## Data frame representation of unmarkedFrame object.
## y.1 y.2 y.3 y.4 trans long type length hab area
## L11 2 1 0 1 L11 50 T 25 OW 2
## L12 0 2 0 0 L12 50 T 25 OW 1
## L124 2 1 1 1 L124 50 T 25 OW 2
## L13 1 2 0 1 L13 50 T 25 OW 1
## L14 14 4 7 2 L14 50 T 25 OW 2
## L15 3 6 5 0 L15 50 S 25 OG 1
## L16 4 1 2 1 L16 50 S 25 OG 2
## L17 8 5 6 4 L17 50 S 25 OG 1
## L95 8 5 2 2 L95 50 S 25 OG 1#Por efecto el análisis fundamental es utilizando ajuste halfnormal
fm1 <- distsamp(~ 1 ~ 1, umf)
summary(fm1) #Revisar la información que se brinda##
## Call:
## distsamp(formula = ~1 ~ 1, data = umf)
##
## Density (log-scale):
## Estimate SE z P(>|z|)
## 5.18 0.134 38.7 0
##
## Detection (log-scale):
## Estimate SE z P(>|z|)
## 2.42 0.132 18.3 1.96e-74
##
## AIC: 172.9927
## Number of sites: 9
## optim convergence code: 0
## optim iterations: 27
## Bootstrap iterations: 0
##
## Survey design: line-transect
## Detection function: halfnorm
## UnitsIn: m
## UnitsOut: hahist(fm1, col=456)
#Ajustando la unidad del resultado a kilómetros cuadrados
(fm1 <- distsamp(~ 1 ~ 1, umf, output="density", unitsOut="kmsq"))##
## Call:
## distsamp(formula = ~1 ~ 1, data = umf, output = "density", unitsOut = "kmsq")
##
## Density:
## Estimate SE z P(>|z|)
## 9.79 0.134 73.2 0
##
## Detection:
## Estimate SE z P(>|z|)
## 2.42 0.132 18.3 1.89e-74
##
## AIC: 172.9927#Ajustando Otros modelos de detección
hn_Null <- distsamp (~1 ~1, umf, keyfun = "halfnorm")
haz_Null <-distsamp (~1 ~1, umf, keyfun = "hazard") #lowest AIC ## [1] "Hessian is singular. Try using fewer covariates or providing starting values."uni_Null <- distsamp (~1 ~1, umf, keyfun = "uniform")
exp_Null <- distsamp (~1 ~1, umf, keyfun = "exp") #Construimos una lista de objetos, para comparar los modelos "AIC"
fmList <- fitList(hn_Null=hn_Null,haz_Null=haz_Null, uni_Null=uni_Null, exp_Null=exp_Null )
modSel(fmList)## Hessian is singular.## nPars AIC delta AICwt cumltvWt
## hn_Null 2 172.99 0.00 1.0e+00 1.00
## exp_Null 2 189.97 16.98 2.1e-04 1.00
## uni_Null 1 189.97 16.98 2.1e-04 1.00
## haz_Null 3 191.97 18.98 7.6e-05 1.00par(mfrow=c(2,2))
hist(hn_Null,col=456, main="halfnorm")
hist(haz_Null, col=456, main="hazard")
hist(uni_Null, col=456, main="uniform")
#Ajustando el "Goodness of fit test" de la detección por covariables
model1 <-distsamp (~1 ~area, umf, keyfun = "halfnorm", unitsOut="ha")
model1##
## Call:
## distsamp(formula = ~1 ~ area, data = umf, keyfun = "halfnorm",
## unitsOut = "ha")
##
## Density:
## Estimate SE z P(>|z|)
## (Intercept) 5.306 0.313 16.976 1.23e-64
## area -0.087 0.198 -0.438 6.61e-01
##
## Detection:
## Estimate SE z P(>|z|)
## 2.42 0.132 18.3 1.89e-74
##
## AIC: 174.7997model2 <-distsamp (~1 ~length, umf, keyfun = "halfnorm",unitsOut="ha")
model2##
## Call:
## distsamp(formula = ~1 ~ length, data = umf, keyfun = "halfnorm",
## unitsOut = "ha")
##
## Density:
## Estimate SE z P(>|z|)
## (Intercept) 0.0681 31598 2.15e-06 1
## length 0.2045 1264 1.62e-04 1
##
## Detection:
## Estimate SE z P(>|z|)
## 2.42 0.132 18.2 2.19e-74
##
## AIC: 174.9927model3 <-distsamp (~1 ~hab+area+length, umf, keyfun = "halfnorm",unitsOut="ha")
model3##
## Call:
## distsamp(formula = ~1 ~ hab + area + length, data = umf, keyfun = "halfnorm",
## unitsOut = "ha")
##
## Density:
## Estimate SE z P(>|z|)
## (Intercept) 0.0857 2.51e+04 3.42e-06 1.00000
## habOW -0.6658 2.14e-01 -3.12e+00 0.00183
## area 0.1506 2.12e-01 7.10e-01 0.47783
## length 0.2079 1.00e+03 2.07e-04 0.99983
##
## Detection:
## Estimate SE z P(>|z|)
## 2.42 0.132 18.3 2.04e-74
##
## AIC: 168.87#Otros métodos para ajustar el modelo y cálculos
backTransform(fm1, type="state") # animals / ha## Backtransformed linear combination(s) of Density estimate(s)
##
## Estimate SE LinComb (Intercept)
## 17784 2379 9.79 1
##
## Transformation: expbackTransform(fm1, type="det") # half-normal SD## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb sigma(Intercept)
## 11.2 1.48 2.42 1
##
## Transformation: exphist(fm1, xlab="Distance (m)") # Only works when there are no det covars
# Empirical Bayes estimates of posterior distribution for N_i
plot(ranef(fm1, K=50))
# Effective strip half-width
(eshw <- integrate(gxhn, 0, 20, sigma=10.9)$value)## [1] 12.7523# Detection probability
eshw / 20 # 20 is strip-width## [1] 0.6376152# Halfnormal. Covariates affecting both density and and detection.
(fm2 <- distsamp(~hab ~ area, umf))##
## Call:
## distsamp(formula = ~hab ~ area, data = umf)
##
## Density:
## Estimate SE z P(>|z|)
## (Intercept) 5.1542 0.323 15.963 2.31e-57
## area 0.0291 0.205 0.142 8.87e-01
##
## Detection:
## Estimate SE z P(>|z|)
## sigma(Intercept) 2.678 0.244 11 5.30e-28
## sigmahabOW -0.486 0.243 -2 4.55e-02
##
## AIC: 171.1749#determining predicted values for each covariate level
m.half.1.hab <- distsamp(~1 ~hab, umf, keyfun="halfnorm", output="density", unitsOut="kmsq")
m.hab <- data.frame(hab=factor(c( "OW", "OG")))
hab.pred <- predict(m.half.1.hab, type="state", newdata=m.hab, appendData=TRUE)
hab.pred ## Predicted SE lower upper hab
## 1 12927.39 2315.776 9099.745 18365.08 OW
## 2 23854.11 3726.867 17562.046 32400.48 OGGraficando mis resultados de las covariables
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.2.3ggplot(hab.pred, aes(x=hab, y=Predicted)) +
geom_bar(position=position_dodge(), stat="identity", colour = 'blue', fill = 'light blue') +
geom_errorbar(aes(ymin=Predicted-SE, ymax=Predicted+SE), width=.2, colour = 'blue') +
scale_x_discrete(name = expression(bold('Vegetation Type'))) +
scale_y_continuous(name = expression(bold(paste("Cantidad / ", km^bold("2"))))) +
theme_bw()
#Tests para la detección del tamaño "Cantidad"
#Plot de la Regresión del cluster size en distance
plot(dists.trunc$dist, dists.trunc$Cantidad)
cluster.det <- lm(dists.trunc$Cantidad ~ dists.trunc$dist)
anova(cluster.det)## Analysis of Variance Table
##
## Response: dists.trunc$Cantidad
## Df Sum Sq Mean Sq F value Pr(>F)
## dists.trunc$dist 1 1.567 1.5675 0.7787 0.3796
## Residuals 103 207.347 2.0131summary(cluster.det) # Si p<0.15, no es menor a este valor, es mejor utilizar una transformación loaritmica para mejorar la incertidumbre ##
## Call:
## lm(formula = dists.trunc$Cantidad ~ dists.trunc$dist)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0891 -0.8913 0.0895 1.0012 5.2093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.65514 0.24042 11.044 <2e-16 ***
## dists.trunc$dist 0.02230 0.02527 0.882 0.38
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.419 on 103 degrees of freedom
## Multiple R-squared: 0.007503, Adjusted R-squared: -0.002133
## F-statistic: 0.7787 on 1 and 103 DF, p-value: 0.3796#Regress log(cluster size) on distance
plot(dists.trunc$dist, log(dists.trunc$Cantidad))
l.herd <- log(dists.trunc$Cantidad)
cluster.det <- lm(l.herd ~ dists.trunc$dist)
summary(cluster.det) ##
## Call:
## lm(formula = l.herd ~ dists.trunc$dist)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0018 -0.2392 0.1595 0.4162 1.1823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.849547 0.090072 9.432 1.36e-15 ***
## dists.trunc$dist 0.007822 0.009468 0.826 0.411
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5315 on 103 degrees of freedom
## Multiple R-squared: 0.006583, Adjusted R-squared: -0.003062
## F-statistic: 0.6825 on 1 and 103 DF, p-value: 0.4106Using this information to convert “Cantidad” density to individual density
bt.est <- exp( 0.849547) #exp of the intercept
bt.est## [1] 2.338587bt.se <- 0.090072/(1/bt.est) #SE of intercept divided by 1 over the back-transformed estimate
bt.se## [1] 0.2106412#Si la regresión del tamaño de las densidades no es fuerte ( p > 0,15 es una regla de oro ) y no hay covariables importantes de la densidad en la parte superior del modelo, entonces el tamaño del grupo medio en general es apropiado para convertir la densidad del grupo a la densidad individuo.
avg.herd <- mean(dists.trunc$Cantidad) #Tamaño esperado del grupo
avg.herd## [1] 2.828571sd.herd <- sd(dists.trunc$Cantidad) #intervalo de confianza del grupo
sd.herd## [1] 1.417318dim(dists.trunc)## [1] 105 4