RSF analysis of mountain goats (Section 4.1) (2024)

S. Muff, J. Signer, J. Fieberg

03 März, 2019

Load libraries and data
glmmTMB()
- Model M1:
- Model M2
- Model M3
- Model M4 (with fixed intercept variance)
- Model M4 (with intercept variance estimated)
INLA (only model M4)
Session Info

Load libraries and data

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.3333333

glmmTMB()

Model M1:

Fixed Effects: elevation
Random Effects: intercept only

goats.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 ' ' 1

Model M2

Fixed Effects: elevation, aspect
Random Effects: intercept only

goats.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 ' ' 1

Model M3

Fixed Effects: elevation, aspect
Random Effects: intercepts and slopes (elevation, aspect)

goats.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 ' ' 1

Model M4 (with fixed intercept variance)

Fixed Effects: elevation
Random Effects: intercept (with large fixed variance), elevation, aspect

Here, 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 ' ' 1

Model M4 (with intercept variance estimated)

For 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 ' ' 1

INLA (only model M4)

Let 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$ID

For 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-11

Since 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.9508333

Source 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.001465

inla_mmarginal(goats.M4.inla)

## Mode of variance for ID2 Mode of variance for ID3 ## 0.1134600 0.7537896

Session Info

devtools::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)