Analysis of mountain goat data from Alaska and Washington

Programmer: John Fieberg

R Code in support of Fieberg et al. Do Capture and Survey Methods Influence Whether Marked Animals are Representative of Unmarked Animals?

# Clear out working memory 
  rm(list = ls())

Mountain goats in Alaska

# set global option for output when producing html file
  opts_chunk$set(comment=NA, fig.width=6, fig.height=6)

load libraries

  library(lme4)
  library(doBy)
  library(gplots)
  library(mosaic)
  library(AICcmodavg)
  library(MuMIn)

Function for getting SE

  SE<-function(x){sqrt(var(x)/length(x))} 

Read in data:

  mtg<-read.csv("../Data/MTG_Sight_Alaska.csv")

Look at number of observations for each animal

  table(mtg$ID)

LG003 LG004 LG005 LG007 LG008 LG009 LG010 LG011 LG015 LG016 LG017 LG020 
    1     1     1     1     1     1     1     1     1     1     1     1 
LG021 LG022 LG023 LG024 LG025 LG026 LG027 LG029 LG030 LG036 LG037 LG038 
    1     1     1     2     1     1     2     2     1     2     2     2 
LG039 LG040 LG041 LG045 LG049 LG050 LG051 LG052 LG053 LG054 LG055 LG056 
    1     3     2     2     1     1     2     2     2     1     1     1 
LG057 LG058 LG059 LG060 LG061 LG062 LG063 LG064 LG066 LG068 LG069 LG070 
    1     1     1     2     3     2     1     2     1     1     2     1 
LG071 LG074 LG075 LG076 LG077 LG078 LG080 LG081 LG083 LG084 LG086 LG087 
    1     1     1     2     2     1     1     1     2     2     2     2 
LG088 LG089 LG090 LG091 LG092 LG093 LG094 LG095 LG096 LG097 LG098 LG099 
    2     2     1     1     1     1     1     1     1     1     1     1 
LG100 LG101 LG102 LG103 LG104 LG105 LG106 LG108 LG109 LG110 LG111 LG112 
    1     1     1     1     1     1     1     1     1     1     1     1 
LG113 LG114 LG115 LG116 LG117 LG118 LG119 LG120 
    1     1     1     1     1     1     1     1 

Distribution of time since capture

  plot(density(mtg$Day.Since.Capture, from=0), xlab="Days since capture", ylab="Density", main="")  

plot of chunk unnamed-chunk-8

Since so few animals were observed more than once, fit a glm (assuming independence), and then use bootstrap for standard errors if necessary (e.g., if statistically significant when independence is assumed).

  fit.mtg<-glm(Seen.~Day.Since.Capture, family=binomial(), data=mtg)
  summary(fit.mtg) 

Call:
glm(formula = Seen. ~ Day.Since.Capture, family = binomial(), 
    data = mtg)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
 -1.19   -1.18   -1.11    1.18    1.25  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)
(Intercept)        0.018938   0.252848    0.07     0.94
Day.Since.Capture -0.000416   0.001044   -0.40     0.69

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 164.89  on 118  degrees of freedom
Residual deviance: 164.73  on 117  degrees of freedom
AIC: 168.7

Number of Fisher Scoring iterations: 3

Could also consider controlling for area if pdetect depeneded on area….

  fit.mtg2<-glm(Seen.~Day.Since.Capture+Area2, family=binomial(), data=mtg)
  summary(fit.mtg2) 

Call:
glm(formula = Seen. ~ Day.Since.Capture + Area2, family = binomial(), 
    data = mtg)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.344  -1.139  -0.961   1.202   1.430  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)
(Intercept)        1.56e+01   1.46e+03    0.01     0.99
Day.Since.Capture -6.21e-04   1.09e-03   -0.57     0.57
Area2East Berners -1.59e+01   1.46e+03   -0.01     0.99
Area2Lions Head   -1.52e+01   1.46e+03   -0.01     0.99
Area2Sinclair     -1.56e+01   1.46e+03   -0.01     0.99
Area2Villard      -1.59e+01   1.46e+03   -0.01     0.99

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 164.89  on 118  degrees of freedom
Residual deviance: 161.41  on 113  degrees of freedom
AIC: 173.4

Number of Fisher Scoring iterations: 14

Days since capture does not appear to be significant. What about “years since capture” (since the number of days since capture is bi-modal)

  mtg$Itime<-I(mtg$Day.Since.Capture > 200) 
  fit.mtg3<-glmer(Seen.~Itime + (1|ID), family=binomial(), nAGQ=10, data=mtg)
  summary(fit.mtg3)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: Seen. ~ Itime + (1 | ID)
   Data: mtg

     AIC      BIC   logLik deviance df.resid 
   170.6    178.9    -82.3    164.6      116 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.013 -1.013 -0.913  0.987  1.095 

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 1.73e-14 1.32e-07
Number of obs: 119, groups:  ID, 92

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.0267     0.2310    0.12     0.91
ItimeTRUE    -0.2090     0.3808   -0.55     0.58

Correlation of Fixed Effects:
          (Intr)
ItimeTRUE -0.607

Test and confidence interval for a difference in proportions

  prop.test(Seen.~Itime, data=mtg)

    2-sample test for equality of proportions with continuity
    correction

data:  t(table_from_formula)
X-squared = 0.129, df = 1, p-value = 0.7195
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.1515  0.2558
sample estimates:
prop 1 prop 2 
0.5455 0.4933 
  with(data=mtg,prop.table(table(Seen., Itime), margin=2))
     Itime
Seen.  FALSE   TRUE
    0 0.4933 0.5455
    1 0.5067 0.4545

Mountain Goats in WA

Read in data:

  mtgwa<-read.csv("../Data/WAgoats.csv")

Fit full model, with interaction between yr.s.cap & method, then model without the interaction

  fit.glmer1<-glmer(observed~grpsize+vegmidpt + terrainob+ method + yr.s.cap + method*yr.s.cap + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer1)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: observed ~ grpsize + vegmidpt + terrainob + method + yr.s.cap +  
    method * yr.s.cap + (1 | id)
   Data: mtgwa

     AIC      BIC   logLik deviance df.resid 
    97.4    119.1    -40.7     81.4      103 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.087  0.016  0.264  0.489  1.051 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.228    0.477   
Number of obs: 111, groups:  id, 52

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)           1.6060     1.2516    1.28    0.199  
grpsize               0.2508     0.1295    1.94    0.053 .
vegmidpt             -0.0113     0.0109   -1.03    0.302  
terrainobY           -0.5268     0.6647   -0.79    0.428  
methodHeli           -0.1484     1.4895   -0.10    0.921  
yr.s.cap             -0.4628     0.7389   -0.63    0.531  
methodHeli:yr.s.cap   0.3297     1.0136    0.32    0.745  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) grpsiz vgmdpt trrnbY mthdHl yr.s.c
grpsize     -0.231                                   
vegmidpt    -0.071  0.090                            
terrainobY  -0.438 -0.132 -0.161                     
methodHeli  -0.601 -0.108 -0.221  0.167              
yr.s.cap    -0.766 -0.034 -0.087  0.131  0.621       
mthdHl:yr..  0.540  0.148  0.145 -0.146 -0.889 -0.743

Drop interaction

  fit.glmer2<-glmer(observed~grpsize+vegmidpt + method + terrainob+ yr.s.cap + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer2)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: observed ~ grpsize + vegmidpt + method + terrainob + yr.s.cap +  
    (1 | id)
   Data: mtgwa

     AIC      BIC   logLik deviance df.resid 
    95.5    114.5    -40.8     81.5      104 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.253  0.017  0.273  0.507  1.105 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.189    0.435   
Number of obs: 111, groups:  id, 52

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   1.3949     1.0233    1.36    0.173  
grpsize       0.2452     0.1265    1.94    0.052 .
vegmidpt     -0.0118     0.0107   -1.10    0.270  
methodHeli    0.2835     0.6767    0.42    0.675  
terrainobY   -0.4988     0.6522   -0.76    0.444  
yr.s.cap     -0.2871     0.4864   -0.59    0.555  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) grpsiz vgmdpt mthdHl trrnbY
grpsize    -0.370                            
vegmidpt   -0.179  0.068                     
methodHeli -0.300  0.050 -0.204              
terrainobY -0.434 -0.119 -0.142  0.089       
yr.s.cap   -0.632  0.111  0.029 -0.163  0.022

Drop other non-significant variables (control only for groupsize) to test for robustness. First, include interaction

  fit.glmer3<-glmer(observed~grpsize+ method + yr.s.cap + method*yr.s.cap+(1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer3)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: observed ~ grpsize + method + yr.s.cap + method * yr.s.cap +  
    (1 | id)
   Data: mtgwa

     AIC      BIC   logLik deviance df.resid 
    95.6    111.8    -41.8     83.6      105 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.189  0.013  0.291  0.521  0.805 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.108    0.329   
Number of obs: 111, groups:  id, 52

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)            1.005      1.057    0.95    0.342  
grpsize                0.261      0.125    2.10    0.036 *
methodHeli            -0.285      1.399   -0.20    0.839  
yr.s.cap              -0.483      0.692   -0.70    0.485  
methodHeli:yr.s.cap    0.398      0.955    0.42    0.677  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) grpsiz mthdHl yr.s.c
grpsize     -0.309                     
methodHeli  -0.641 -0.078              
yr.s.cap    -0.815 -0.032  0.608       
mthdHl:yr..  0.552  0.141 -0.888 -0.726

Now, drop interaction

  fit.glmer3b<-glmer(observed~grpsize+ method + yr.s.cap +(1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer3b)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: observed ~ grpsize + method + yr.s.cap + (1 | id)
   Data: mtgwa

     AIC      BIC   logLik deviance df.resid 
    93.7    107.3    -41.9     83.7      106 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.433  0.015  0.287  0.536  0.741 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.095    0.308   
Number of obs: 111, groups:  id, 52

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)    0.768      0.854    0.90    0.368  
grpsize        0.255      0.123    2.08    0.038 *
methodHeli     0.236      0.641    0.37    0.713  
yr.s.cap      -0.277      0.469   -0.59    0.554  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) grpsiz mthdHl
grpsize    -0.477              
methodHeli -0.372  0.106       
yr.s.cap   -0.709  0.104 -0.157

Now, method and yr.s.cap as only predictors (along with groupsize)

  fit.glmer4<-glmer(observed~grpsize+ method + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer4)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: observed ~ grpsize + method + (1 | id)
   Data: mtgwa

     AIC      BIC   logLik deviance df.resid 
    92.1    102.9    -42.1     84.1      107 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1374  0.0131  0.2939  0.5505  0.6739 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.126    0.355   
Number of obs: 111, groups:  id, 52

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)    0.413      0.606    0.68    0.496  
grpsize        0.265      0.125    2.12    0.034 *
methodHeli     0.183      0.638    0.29    0.774  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) grpsiz
grpsize    -0.573       
methodHeli -0.690  0.126
  fit.glmer5<-glmer(observed~grpsize+  yr.s.cap + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer5)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: observed ~ grpsize + yr.s.cap + (1 | id)
   Data: mtgwa

     AIC      BIC   logLik deviance df.resid 
    91.9    102.7    -41.9     83.9      107 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.280  0.015  0.302  0.561  0.729 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.0255   0.16    
Number of obs: 111, groups:  id, 52

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)    0.887      0.787    1.13    0.260  
grpsize        0.251      0.121    2.08    0.037 *
yr.s.cap      -0.250      0.459   -0.54    0.586  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
         (Intr) grpsiz
grpsize  -0.465       
yr.s.cap -0.840  0.122

Yr.s.cap is not even close to significant in any of the models. Construct plots to look at effect size.

AK mountain goats: basic proportions and confidence intervals for proportion seen vs. time since cap

  mtg.ak<-summaryBy(Seen.~Itime, FUN=c(mean, SE),data=mtg) 

Model based predictions for Alaskan mtgs

  xmat<-matrix(rbind(c(1,0), c(1,1)), 2,2)
  mt.ak.prd<- matrix(NA,2,3)
  mt.ak.prd[,1]<-plogis(xmat%*%fixef(fit.mtg3))
  se.prd<-diag(xmat%*%vcov(fit.mtg3)%*% t(xmat)) 
  mt.ak.prd[,3]<- plogis(xmat%*%fixef(fit.mtg3)+1.96*se.prd)
  mt.ak.prd[,2]<- plogis(xmat%*%fixef(fit.mtg3)-1.96*se.prd)

Now for WA mtgs. Use fit.glmer5, with grpsize and yr.s.cap (lowest AIC of models containing yr.s.cap) Get predictions for median group size (which varied a bit by yr.s.cap, but not that much)

  tapply(mtgwa$grpsize, mtgwa$yr.s.cap, mean)
     0      1      2      3 
 7.000  9.065  8.821 46.000 
  tapply(mtgwa$grpsize, mtgwa$yr.s.cap, median)
 0  1  2  3 
 6  4  3 46 
  mgs<-median(mtgwa$grpsize)
  xmat<-matrix(rbind(c(1,mgs,0), c(1,mgs,1), c(1,mgs,2)), 3,3)
  mt.wa.prd<- matrix(NA,3,3)
  mt.wa.prd[,1]<-plogis(xmat%*%fixef(fit.glmer5))
  se.prd.wa<-diag(xmat%*%vcov(fit.glmer5)%*% t(xmat)) 
  mt.wa.prd[,3]<- plogis(xmat%*%fixef(fit.glmer5)+1.96*se.prd.wa)
  mt.wa.prd[,2]<- plogis(xmat%*%fixef(fit.glmer5)-1.96*se.prd.wa)

Now, get proportion observed as function of time since cap

  mtgwa$Iobs<-I(mtgwa$observed=="Y")
  mtg.wa.ob<-summaryBy(Iobs~yr.s.cap,FUN=c(mean, SE), data=mtgwa)

# Now, plot with helicopter/ground as well for washington data, using glmer3 (with interaction) between method & yr since cap
  xmat<-matrix(rbind(c(1,mgs,0,0,0),c(1,mgs,0,1,0),c(1,mgs,0,2,0), 
                    c(1,mgs,1,0,1),c(1,mgs,1,1,1),c(1,mgs,1,2,1)), 6,5)
  mt.wa.prd<- matrix(NA,6,3)
  mt.wa.prd[,1]<-plogis(xmat%*%fixef(fit.glmer3))
  se.prd.wa<-diag(xmat%*%vcov(fit.glmer3)%*% t(xmat)) 
  mt.wa.prd[,3]<- plogis(xmat%*%fixef(fit.glmer3)+1.96*se.prd.wa)
  mt.wa.prd[,2]<- plogis(xmat%*%fixef(fit.glmer3)-1.96*se.prd.wa)

Plot results

par(mfrow=c(1,1), bty="L", cex.lab=1.8, cex=1, cex.axis=1.8, mar=c(6,6,4.1,2.1) , xpd=T, oma=c(0,1,0,0)) 
  plotCI(c(1,2), mt.ak.prd[,1], ui=mt.ak.prd[,3], li=mt.ak.prd[,2], xlab="Years since capture", ylab="Probability of detection", ylim=c(0,1), xlim=c(0,2), gap=0, pch=17, cex=2,xaxt="n")
  axis(side=1, c(0,1,2))
  plotCI(c(c(0,1,2)-.02, c(0,1,2)+.02), mt.wa.prd[,1], ui=mt.wa.prd[,3], li=mt.wa.prd[,2], xlab="", ylab="",add=T, ylim=c(0,1), xlim=c(0,2), gap=0, pch=c(rep(16,3), rep(1,3)),cex=2, xaxt="n")   
  legend("bottomright", pch=c(17,16,1), cex=2,lty=1, title="", legend=c("Alaska Helicopter/Helicopter", "Washington Ground/Helicopter", "Washington Helicopter/Helicopter"), bty="n")

plot of chunk unnamed-chunk-23

For WA, where there were no missing values, also calculate model-averaged coefficients to get at effect sizes (use all models except those with interactions between method and years since capture; these models had larger AIC's anyway)

  Cands <- list(fit.glmer2, fit.glmer3b, fit.glmer4, fit.glmer5)

model selection

  aictab(cand.set = Cands)
Warning: 
Model names have been supplied automatically in the table

Model selection based on AICc :

     K  AICc Delta_AICc AICcWt Cum.Wt     LL
Mod4 4 92.27       0.00   0.42   0.42 -41.94
Mod3 4 92.48       0.21   0.38   0.80 -42.05
Mod2 5 94.32       2.06   0.15   0.95 -41.87
Mod1 7 96.59       4.32   0.05   1.00 -40.75

model average for difference between period 1 vs period 4

  summary(mod.avg(fit.glmer2, fit.glmer3b, fit.glmer4, fit.glmer5))
Error: error in evaluating the argument 'object' in selecting a method for function 'summary': Error: could not find function "mod.avg"

Effect sizes:

  c(exp(-0.161),  exp(-0.161-1.96*0.46395), exp(-0.161+1.96*0.46395))
[1] 0.8513 0.3429 2.1135
  c(exp(0.1188361),exp(0.1188361-1.96*0.4996),exp(0.1188361+1.96*0.4996))
[1] 1.126 0.423 2.998
# Print out Session info
  sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] RMark_2.1.8     MuMIn_1.10.5    AICcmodavg_2.00 gplots_2.14.1  
 [5] knitr_1.6       mosaic_0.9.1-3  lattice_0.20-29 dplyr_0.2      
 [9] car_2.0-21      doBy_4.5-10     MASS_7.3-33     survival_2.37-7
[13] nlme_3.1-117    ggplot2_1.0.0   lme4_1.1-7      Rcpp_0.11.2    
[17] Matrix_1.1-4   

loaded via a namespace (and not attached):
 [1] assertthat_0.1     bitops_1.0-6       caTools_1.17      
 [4] coda_0.16-1        colorspace_1.2-4   digest_0.6.4      
 [7] evaluate_0.5.5     expm_0.99-1.1      formatR_1.0       
[10] gdata_2.13.3       ggdendro_0.1-14    grid_3.1.1        
[13] gridExtra_0.9.1    gtable_0.1.2       gtools_3.4.1      
[16] htmltools_0.2.4    KernSmooth_2.23-12 labeling_0.3      
[19] markdown_0.7.4     matrixcalc_1.0-3   mime_0.1.2        
[22] minqa_1.2.3        mosaicData_0.9.1   msm_1.4           
[25] munsell_0.4.2      mvtnorm_1.0-0      nloptr_1.0.4      
[28] nnet_7.3-8         parallel_3.1.1     plyr_1.8.1        
[31] proto_0.3-10       reshape2_1.4       rmarkdown_0.2.49  
[34] scales_0.2.4       snowfall_1.84-6    stringr_0.6.2     
[37] tools_3.1.1        yaml_2.1.13