#Generalized Linear Mixed Effects Regression
###Before you begin
-Install the lme4 package
-Download the dataset for this exercise and save it to your working directory: final_winter_modeldata2.csv
###Read the dataset into R
winter_deer_data1=read.csv('final_winter_modeldata2.csv') # read data file into RThese data were obtained by extracting values for each GPS collar location from GIS layers for the area these deer were occupying. Random points were created and had values extracted to classify what was “available” to the deer. The data has a “used” column that is coded 1 or 0 for used or unused respectively. Deer ID is “altid”.
###Convert catagorical variables to factors
Catgorical variables must be converted to factors to be recognized properly in R for this exercise
winter_deer_data1$altid=factor(winter_deer_data1$altid) #convert catagorical variable to a vector of factor variables
winter_deer_data1$veg_class=factor(winter_deer_data1$veg_class) #convert catagorical variable to a vector of factor variables###Standardize all continuous variables
To avoid convergence issues related to scaling, it is important to standardize continuous variables
stand_dist_to_water=(winter_deer_data1$dist_to_water-mean(winter_deer_data1$dist_to_water))/sd(winter_deer_data1$dist_to_water) #standardize continuous variable
stand_cos_aspect=(winter_deer_data1$cos_aspect-mean(winter_deer_data1$cos_aspect))/sd(winter_deer_data1$cos_aspect) #standardize continuous variable
stand_sin_aspect=(winter_deer_data1$sin_aspect-mean(winter_deer_data1$sin_aspect))/sd(winter_deer_data1$sin_aspect) #standardize continuous variable
stand_elevation=(winter_deer_data1$elevation-mean(winter_deer_data1$elevation))/sd(winter_deer_data1$elevation) #standardize continuous variable
stand_slope=(winter_deer_data1$slope-mean(winter_deer_data1$slope))/sd(winter_deer_data1$slope) #standardize continuous variable####Load the lme4 package
library(lme4) #load lme4 package## Loading required package: Matrix##Build a general mixed-effect additive model with random intercept
Here we build a basic generalized linear mixed-effects model. We use the glmer function and reference the used column from our data to compare used vs available points. This is an additive model that accounts for random effects allowing for random slope by individual. Note: “altid” is the identification number for each deer. The values associated with each GPS location have the deer ID value. This is treated as the “group” for mixed-effects in this example.
# generalized linear mixed effect additive model with random intercept for individual
winter_pequop_model=glmer(used~stand_dist_to_water+stand_cos_aspect
+stand_elevation+stand_slope+veg_class+
(1|altid),family="binomial",data=winter_deer_data1,
nAGQ = 1, na.action="na.fail")
summary(winter_pequop_model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## used ~ stand_dist_to_water + stand_cos_aspect + stand_elevation +
## stand_slope + veg_class + (1 | altid)
## Data: winter_deer_data1
##
## AIC BIC logLik deviance df.resid
## 14993.5 15081.5 -7485.8 14971.5 22061
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -27.6196 -0.3260 0.0032 0.3248 7.6500
##
## Random effects:
## Groups Name Variance Std.Dev.
## altid (Intercept) 0.9835 0.9917
## Number of obs: 22072, groups: altid, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.41433 0.14108 2.94 0.00332 **
## stand_dist_to_water -0.74022 0.02806 -26.38 < 2e-16 ***
## stand_cos_aspect -0.18377 0.02122 -8.66 < 2e-16 ***
## stand_elevation 1.56434 0.04095 38.20 < 2e-16 ***
## stand_slope 1.22337 0.03380 36.20 < 2e-16 ***
## veg_class3 -0.57140 0.05883 -9.71 < 2e-16 ***
## veg_class4 -4.67309 0.76972 -6.07 1.27e-09 ***
## veg_class6 -0.05232 0.12770 -0.41 0.68199
## veg_class7 1.71760 0.13837 12.41 < 2e-16 ***
## veg_class8 -2.67071 0.45203 -5.91 3.46e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) stn___ stnd__ stnd_l stnd_s vg_cl3 vg_cl4 vg_cl6 vg_cl7
## stnd_dst_t_ 0.034
## stnd_cs_spc -0.009 -0.002
## stand_elvtn 0.075 0.207 0.009
## stand_slope 0.033 -0.219 0.057 -0.420
## veg_class3 -0.146 -0.054 0.051 -0.197 -0.241
## veg_class4 -0.016 0.033 0.001 -0.036 -0.026 0.037
## veg_class6 -0.073 0.038 0.025 -0.284 0.058 0.255 0.024
## veg_class7 -0.054 -0.169 0.008 0.090 -0.031 0.158 0.000 0.048
## veg_class8 -0.032 -0.006 -0.028 -0.149 0.007 0.104 0.010 0.097 0.017The summary will give the coefficients for each parameter evaluated in the model. Positive coefficients mean selection for that variable with the exception of “distance to” variables" which are opposite. The summary also provides an AIC value to be used for comparison to models that will be built in the below excercises. Keep an eye on the standard error and significance values as you build models for anything that looks out of the ordinary.
***Note: in the following exercises it is important to understand that convergence issues may occur as models become too complex or if there is large variation in your group effect. Think about the relationships that may exist between covariates and omit those that may be less important after evaluating significance from several models.
##Exercise 1
###Building more complex models
Explore building more complex models. Generalized linear mixed effect models are flexible allowing the user to specify random slope, random intercept and other random effects. Be wary in blindly building models without biological reasoning. Below is a list of coding options for random effects code modifications:
(1|group) = random group intercept
(x|group) = random slope of x within group with correlated intercept
(0+x|group) = random slope of x within group: no variation in intercept
(1|group) + (0+x|group) = uncorrelated random intercept and random slope within group
(1|site/block) = intercept varying among sites and among blocks within sites (nested random effects)
site+(1|site:block) = fixed effect of sites plus random variation in intercept among blocks within sites
(x|site/block) = slope and intercept varying among sites and among blocks within sites
(x*site+(x|site:block)) = fixed effect variation of slope and intercept varying among sites and random variation of slope and intercept among blocks within sites
(1|group1)+(1|group2) = intercept varying among crossed random effects (e.g. site, year)
####Example of a model allowing coefficient for elevation to vary among each altid (deer)
winter_pequop_model2=glmer(used~stand_dist_to_water+stand_cos_aspect
+(stand_elevation|altid)+stand_slope+veg_class,
family="binomial",data=winter_deer_data1,
nAGQ = 1, na.action="na.fail")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00376095 (tol =
## 0.001, component 1)summary(winter_pequop_model2) ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## used ~ stand_dist_to_water + stand_cos_aspect + (stand_elevation |
## altid) + stand_slope + veg_class
## Data: winter_deer_data1
##
## AIC BIC logLik deviance df.resid
## 12744.3 12840.4 -6360.2 12720.3 22060
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -264.989 -0.237 0.000 0.248 4.026
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## altid (Intercept) 75.44 8.686
## stand_elevation 116.54 10.795 0.98
## Number of obs: 22072, groups: altid, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.49577 0.30814 -4.854 1.21e-06 ***
## stand_dist_to_water -0.64301 0.02993 -21.482 < 2e-16 ***
## stand_cos_aspect -0.08085 0.02214 -3.652 0.000261 ***
## stand_slope 0.87497 0.03545 24.681 < 2e-16 ***
## veg_class3 -0.82505 0.06128 -13.464 < 2e-16 ***
## veg_class4 -5.33809 0.93785 -5.692 1.26e-08 ***
## veg_class6 0.23839 0.12385 1.925 0.054251 .
## veg_class7 1.34350 0.16274 8.256 < 2e-16 ***
## veg_class8 -1.60572 0.57540 -2.791 0.005261 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) stn___ stnd__ stnd_s vg_cl3 vg_cl4 vg_cl6 vg_cl7
## stnd_dst_t_ -0.010
## stnd_cs_spc -0.002 0.006
## stand_slope 0.056 -0.224 0.021
## veg_class3 -0.059 -0.067 0.036 -0.183
## veg_class4 -0.003 0.027 -0.001 -0.007 0.038
## veg_class6 -0.013 0.028 0.034 0.059 0.265 0.020
## veg_class7 -0.039 -0.149 0.000 -0.010 0.172 0.011 0.061
## veg_class8 0.003 0.002 -0.008 -0.002 0.077 0.004 0.094 0.003
## convergence code: 0
## Model failed to converge with max|grad| = 0.00376095 (tol = 0.001, component 1)##Exercise 2
###Explore interactions between covariates
Interactions within the models are specified using the * symbol. Think about biological interactions that may occur between the variables that have been measured.
##Exercise 3
###Explore quadratic relationships
Try exploring quadratic relationships for covariates. Covariates may not always have a linear relationship with resource selection.
***Note: a quadratic relationship can be indicated by I(covariate^2)