library(ResourceSelection)
library(glmmTMB)
library(INLA)
data(goats)
str(goats)## 'data.frame': 19014 obs. of 8 variables:
## $ STATUS : int 1 1 1 1 1 1 1 1 1 1 ...
## $ ID : int 1 1 1 1 1 1 1 1 1 1 ...
## $ ELEVATION: int 651 660 316 334 454 343 429 493 400 442 ...
## $ SLOPE : num 38.5 39.7 20.5 34.1 41.6 ...
## $ ET : num 35.4 70.7 50 35.4 25 ...
## $ ASPECT : num 243 270 279 266 258 ...
## $ HLI : num 0.918 0.884 0.713 0.864 0.935 ...
## $ TASP : num 0.947 0.699 0.575 0.745 0.829 ...Scale and center variables
goats$ELEVATION <- scale(goats$ELEVATION)
goats$TASP <- scale(goats$TASP)Use and available data by animal
with(goats, prop.table(table(ID, STATUS), 1))## STATUS
## ID 0 1
## 1 0.6666667 0.3333333
## 2 0.6666667 0.3333333
## 3 0.6666667 0.3333333
## 4 0.6666667 0.3333333
## 5 0.6666667 0.3333333
## 6 0.6666667 0.3333333
## 7 0.6666667 0.3333333
## 8 0.6666667 0.3333333
## 9 0.6666667 0.3333333
## 10 0.6666667 0.3333333goats.M1 <- glmmTMB(STATUS ~ ELEVATION + (1|ID), family=binomial(), data = goats)
summary(goats.M1)## Family: binomial ( logit )
## Formula: STATUS ~ ELEVATION + (1 | ID)
## Data: goats
##
## AIC BIC logLik deviance df.resid
## 24199.8 24223.4 -12096.9 24193.8 19011
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.007625 0.08732
## Number of obs: 19014, groups: ID, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.68924 0.03192 -21.592 <2e-16 ***
## ELEVATION 0.11844 0.04634 2.556 0.0106 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1goats.M2 <- glmmTMB(STATUS ~ ELEVATION + TASP + (1|ID), family=binomial(), data = goats)
summary(goats.M2)## Family: binomial ( logit )
## Formula: STATUS ~ ELEVATION + TASP + (1 | ID)
## Data: goats
##
## AIC BIC logLik deviance df.resid
## 23331.1 23362.6 -11661.6 23323.1 19010
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.01303 0.1142
## Number of obs: 19014, groups: ID, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.74221 0.03970 -18.694 < 2e-16 ***
## ELEVATION 0.14362 0.02579 5.568 2.58e-08 ***
## TASP 0.51913 0.01883 27.566 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1goats.M3 <- glmmTMB(STATUS ~ TASP + ELEVATION + (1|ID) + (0+ELEVATION |ID) + (0+TASP|ID), family=binomial(), data = goats)
summary(goats.M3)## Family: binomial ( logit )
## Formula:
## STATUS ~ TASP + ELEVATION + (1 | ID) + (0 + ELEVATION | ID) +
## (0 + TASP | ID)
## Data: goats
##
## AIC BIC logLik deviance df.resid
## 21978.5 22025.7 -10983.3 21966.5 19008
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.9572 0.9783
## ID.1 ELEVATION 1.4011 1.1837
## ID.2 TASP 0.1043 0.3229
## Number of obs: 19014, groups: ID, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.39285 0.31158 -1.261 0.207
## TASP 0.66416 0.10565 6.286 3.25e-10 ***
## ELEVATION 0.06934 0.37622 0.184 0.854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here, we also use a weighted likelihood. To this end, we need to create a variable for the weights, where used points (STATUS=1) keep weight 1, and available points (STATUS=0) obtain a large weight \(W\) (here \(W=1000\)):
goats$weight <- 1000^(1-goats$STATUS)We fit the same model as under M3, again using glmmTMB. Note that we have to manually fix the variance of the intercept first. Start by setting up the model, but do not yet fit it:
goats.M4.tmp <- glmmTMB(STATUS ~ ELEVATION + TASP + (1|ID) + (0+ELEVATION |ID) + (0+TASP|ID), family=binomial(), data = goats,doFit=F, weights = weight)Then fix the standard deviation of the first random term, which is the (1|ID) component in the above model equation. We use \(\sigma=10^3\), which corresponds to a variance of \(10^6\):
goats.M4.tmp$parameters$theta[1] = log(1e3)We need to tell glmmTMB not to change the first entry of the vector of variances, and give all other variances another indicator to make sure they can be freely estimated:
goats.M4.tmp$mapArg = list(theta=factor(c(NA,1:2)))Then fit the model and look at the results:
goats.M4 <- glmmTMB:::fitTMB(goats.M4.tmp)
summary(goats.M4)## Family: binomial ( logit )
## Formula:
## STATUS ~ ELEVATION + TASP + (1 | ID) + (0 + ELEVATION | ID) +
## (0 + TASP | ID)
## Data: goats
## Weights: weight
##
## AIC BIC logLik deviance df.resid
## 105999.6 106038.9 -52994.8 105989.6 19009
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.000e+06 1000.0000
## ID.1 ELEVATION 9.257e-01 0.9621
## ID.2 TASP 1.221e-01 0.3495
## Number of obs: 19014, groups: ID, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.633e-05 3.163e+02 0.000 1.000
## ELEVATION 1.172e-01 3.055e-01 0.384 0.701
## TASP 6.505e-01 1.127e-01 5.770 7.92e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For comparison, we again fit model M4, but without fixing the intercept variance, letting it be estimated instead. Importantly, estimating the intercept variance is the current standard procedure. For this particular RSF case, it does not lead to a real difference, as expected due to the many observations per individual. This confirms that the decision to fix or estimate the intercept variance is not critical for RSFs, in contrast to SSFs (see Discussion in the paper).
goats.M4.2 <- glmmTMB(STATUS ~ ELEVATION + TASP + (1|ID) + (0+ELEVATION |ID) + (0+TASP|ID), family=binomial(), data = goats, weights = weight)
summary(goats.M4.2)## Family: binomial ( logit )
## Formula:
## STATUS ~ ELEVATION + TASP + (1 | ID) + (0 + ELEVATION | ID) +
## (0 + TASP | ID)
## Data: goats
## Weights: weight
##
## AIC BIC logLik deviance df.resid
## 105869.2 105916.3 -52928.6 105857.2 19008
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.6448 0.8030
## ID.1 ELEVATION 0.9151 0.9566
## ID.2 TASP 0.1208 0.3476
## Number of obs: 19014, groups: ID, 10
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.3852 0.2554 -28.918 < 2e-16 ***
## ELEVATION 0.1177 0.3037 0.387 0.698
## TASP 0.6501 0.1121 5.797 6.73e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Let us now carry the analysis of model M4 with random intercept \(\mathsf{N}(0,\sigma_{ID}^2)\) and fixed variance \(\sigma_{ID}^2=10^6\) using INLA. A peculiarity of INLA is that the same variable cannot be used more than once. So for ID we need to generate two new (but identical) variables:
goats$ID2 <- goats$ID3 <- goats$IDFor the fixed effects we use the INLA (default) priors \(\beta \sim \mathsf{N}(0,\sigma_\beta^2)\) with \(\sigma_\beta^2=10^4\). The precisions of the priors are thus set to:
prec.beta.TASP <- 1e-4
prec.beta.ELEVATION <- 1e-4 The INLA formula with the fixed effects TASP and ELEVATION, plus three random effects: one for the individual-specific intercept, and two random slopes for TASP and ELEVATION. Note that the precision (thus \(1/\sigma^2\)) for ID is fixed (fixed=T) at the value of \(10^{-6}\) (thus the variance is fixed at \(10^6\)). The precisions for the random slopes for TASP and ELEVATION are given PC(1,0.05) priors:
formula.inla <-STATUS ~ TASP + ELEVATION +
f(ID,model="iid",hyper=list(theta = list(initial=log(1e-6),fixed=T))) +
f(ID2,TASP,values=1:10,model="iid",
hyper=list(theta=list(initial=log(1),fixed=F,prior="pc.prec",param=c(1,0.05)))) +
f(ID3,ELEVATION,values=1:10,model="iid",
hyper=list(theta=list(initial=log(1),fixed=F,prior="pc.prec",param=c(1,0.05)))) The actual INLA call is then given as follows:
inla.setOption(enable.inla.argument.weights=TRUE)
goats.M4.inla <- inla(formula.inla, family ="binomial", data=goats, weights=goats$weight,
control.fixed = list(
mean = 0,
prec = list(TASP = prec.beta.TASP,
ELEVATION = prec.beta.ELEVATION)
)
)The summary for the posterior distribution of the fixed effects is given as follows:
goats.M4.inla$summary.fixed## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -7.6121409 316.0132621 -628.0527313 -7.6210401 612.3106604
## TASP 0.6512586 0.1319021 0.3898437 0.6502860 0.9177524
## ELEVATION 0.1171479 0.3187987 -0.5182517 0.1170705 0.7523008
## mode kld
## (Intercept) -7.6121409 0.000000e+00
## TASP 0.6483680 2.452803e-12
## ELEVATION 0.1169485 1.996053e-11Since variances are parameterized and treated as precisions, the summary of the respective posterior distributions is given for the precisions:
goats.M4.inla$summary.hyperpar## mean sd 0.025quant 0.5quant 0.975quant
## Precision for ID2 7.577611 3.7246041 2.5240783 6.875620 16.824893
## Precision for ID3 1.177238 0.4822806 0.4720444 1.101542 2.333437
## mode
## Precision for ID2 5.5230223
## Precision for ID3 0.9508333Source R functions for calculating posterior means and medians of the precisions.
source("inla_emarginal.R")
source("inla_mmarginal.R")
inla_emarginal(goats.M4.inla)## Mean of variance for ID2 Mean of variance for ID3
## 0.166317 1.001465inla_mmarginal(goats.M4.inla)## Mode of variance for ID2 Mode of variance for ID3
## 0.1134600 0.7537896devtools::session_info()## Session info -------------------------------------------------------------## setting value
## version R version 3.4.4 (2018-03-15)
## system x86_64, linux-gnu
## ui X11
## language de_CH:de
## collate de_CH.UTF-8
## tz Europe/Zurich
## date 2019-03-03## Packages -----------------------------------------------------------------## package * version date source
## backports 1.1.2 2017-12-13 cran (@1.1.2)
## base * 3.4.4 2018-03-16 local
## compiler 3.4.4 2018-03-16 local
## datasets * 3.4.4 2018-03-16 local
## Deriv 3.8.5 2018-06-11 CRAN (R 3.4.4)
## devtools 1.13.4 2017-11-09 CRAN (R 3.4.2)
## digest 0.6.18 2018-10-10 cran (@0.6.18)
## evaluate 0.10.1 2017-06-24 CRAN (R 3.4.1)
## glmmTMB * 0.2.0 2017-12-11 CRAN (R 3.4.3)
## graphics * 3.4.4 2018-03-16 local
## grDevices * 3.4.4 2018-03-16 local
## grid 3.4.4 2018-03-16 local
## htmltools 0.3.6 2017-04-28 CRAN (R 3.4.0)
## INLA * 18.07.11 2018-07-12 local
## knitr 1.17 2017-08-10 CRAN (R 3.4.1)
## lattice 0.20-35 2017-03-25 CRAN (R 3.4.0)
## lme4 1.1-20 2019-02-04 cran (@1.1-20)
## magrittr 1.5 2014-11-22 CRAN (R 3.4.0)
## MASS 7.3-47 2017-04-21 CRAN (R 3.4.2)
## Matrix * 1.2-14 2018-04-09 CRAN (R 3.4.4)
## memoise 1.1.0 2017-04-21 CRAN (R 3.4.2)
## methods * 3.4.4 2018-03-16 local
## minqa 1.2.4 2014-10-09 CRAN (R 3.4.0)
## nlme 3.1-137 2018-04-07 CRAN (R 3.4.4)
## nloptr 1.2.1 2018-10-03 cran (@1.2.1)
## Rcpp 1.0.0 2018-11-07 cran (@1.0.0)
## ResourceSelection * 0.3-2 2017-02-28 CRAN (R 3.4.1)
## rmarkdown 1.7 2017-11-10 CRAN (R 3.4.2)
## rprojroot 1.3-2 2018-01-03 cran (@1.3-2)
## sp * 1.2-5 2017-06-29 cran (@1.2-5)
## splines 3.4.4 2018-03-16 local
## stats * 3.4.4 2018-03-16 local
## stringi 1.2.4 2018-07-20 cran (@1.2.4)
## stringr 1.3.1 2018-05-10 cran (@1.3.1)
## TMB 1.7.14 2018-06-23 CRAN (R 3.4.4)
## tools 3.4.4 2018-03-16 local
## utils * 3.4.4 2018-03-16 local
## withr 2.1.2 2018-03-15 CRAN (R 3.4.4)
## yaml 2.1.14 2016-11-12 CRAN (R 3.4.1)