library(tidyverse)
library(performance) #for model testing
library(broom) #contains augment
library(car) #to check collinearity
library(patchwork) #multiple graphs on same output
library(ggeffects) #to get modeled means out for plotting
library(ggridges) #for cool histograms
library(Matrix)
library(lme4)
library(datawizard)
library(data.table)
library(broom.mixed)
library(scales)SDS300CR lab 4: mixed models and model selection
Some really helpful tutorials for understanding this lab my tutorial Our Coding Club
1.) Get the data. We will use coral recovery/resistance data from my 2021 publication from github for this lab.
There are A LOT of variables here. We are interested in calculated.recovery.rate at the response variable. This is a rate in % coral cover per year (%/yr) that a reef recovers after a disturbance event. Disturbance events are acute events (like a storm) that cause the coral cover to drop. A reef is considered more resilient if it recovers coral cover faster. So higher % cover is good, lower is not as good.
We are interested in the following explanatory variables:
1. region (what ocean basin is the coral from?)
2. disturbance (type of disturbance that occurred)
3. pre.dist.cc (the coral cover at the site before the disturbance in %)
4. recovery_time (time in years over which recovery rate is measured)
5. hii100km (the human influence index in a 100km radius of the site)
6. distance_to_shore_m (the distance to the nearest major shoreline in meters)
7. dist_to_riv_m (the distance to the nearest major river mount in meters)
8. grav_NC (market gravity - proximity to fish market. This is actually the number of people within a 5km range of the reef)
9. cml_scr (cumulative WCS score– a cumulative measure of human influence - similar to HII, but updated to be reef specific)
10. study (the published paper that these data came from)
a.) Read in the data, check shape and size, and keep only the columns listed above. If YOU are interested in another column or group of columns you can keep them as well.
#read data in here. Check df shape and size.
coral_recovery_data <-read.csv("https://raw.githubusercontent.com/jbaumann3/Baumann-et-al-coral-recovery-/main/data/coral_recov_master_2021.csv")
coral_recovery_data <- coral_recovery_data %>% select("region", "disturbance", "pre.dist.cc", "recovery_time" , "hii100km", "distance_to_shore_m", "dist_to_riv_m", "grav_NC", "cml_scr", "study", "calculated.recovery.rate")
head(coral_recovery_data) region disturbance pre.dist.cc recovery_time hii100km
1 Caribbean Bleaching 32.550 1 259587
2 Caribbean Disease, Dredging 30.000 1 259587
3 Caribbean Bleaching 31.500 NA 123407
4 W. Pacific Storm 35.676 NA 42516
5 Caribbean Bleaching 15.600 2 26030
6 Caribbean Bleaching 20.200 2 100868
distance_to_shore_m dist_to_riv_m grav_NC cml_scr
1 480 310479.8 0.2909292 0.4774284
2 480 310479.8 0.2909292 0.4774284
3 2630 616852.4 0.4806241 0.4352647
4 0 1086527.1 0.4636924 0.4059301
5 2630 233058.3 0.4514431 0.5589611
6 532 201198.5 0.5063026 0.6129730
study calculated.recovery.rate
1 Arias-Gonzalez et al (2017) 4.32
2 Arias-Gonzalez et al (2017) 4.81
3 Manzello (2018) NA
4 Madin et al (2018) NA
5 Buglass et al (2015) -2.10
6 Buglass et al (2015) 0.55dim(coral_recovery_data)[1] 223 11b.) For each of our numerical explanatory variables, we need to center and scale them so that we standardize them. This is a good practice for complex model building! Center will make the data normal and scale will change the distribution to make sd=1. Use these scaled data for your models!
#example: df$colname<-scale(df$colname, center=TRUE, scale=TRUE)
#You CAN make a new column if you want to retain the original unscaled data!
coral_recovery_data$pre.dist.cc <- scale(coral_recovery_data$pre.dist.cc, center=TRUE, scale=TRUE)
coral_recovery_data$hii100km <- scale(coral_recovery_data$hii100km, center=TRUE, scale=TRUE)
coral_recovery_data$recovery_time <- scale(coral_recovery_data$recovery_time, center=TRUE, scale=TRUE)
coral_recovery_data$distance_to_shore_m <- scale(coral_recovery_data$distance_to_shore_m, center=TRUE, scale=TRUE)
coral_recovery_data$dist_to_riv_m <- scale(coral_recovery_data$dist_to_riv_m, center=TRUE, scale=TRUE)
coral_recovery_data$grav_NC <- scale(coral_recovery_data$grav_NC, center=TRUE, scale=TRUE)
coral_recovery_data$cml_scr <- scale(coral_recovery_data$cml_scr, center=TRUE, scale=TRUE)
head(coral_recovery_data) region disturbance pre.dist.cc recovery_time hii100km
1 Caribbean Bleaching -0.22330608 -1.0663112 1.7127547
2 Caribbean Disease, Dredging -0.36202416 -1.0663112 1.7127547
3 Caribbean Bleaching -0.28042528 NA 0.4115162
4 W. Pacific Storm -0.05325403 NA -0.3614202
5 Caribbean Bleaching -1.14537331 -0.7993214 -0.5189487
6 Caribbean Bleaching -0.89513678 -0.7993214 0.1961497
distance_to_shore_m dist_to_riv_m grav_NC cml_scr
1 -0.2723830 -0.8208457 -1.3528549 -0.04761585
2 -0.2723830 -0.8208457 -1.3528549 -0.04761585
3 -0.2097209 -0.5676273 -0.6143039 -0.23516921
4 -0.2863727 -0.1794390 -0.6802252 -0.36565605
5 -0.2097209 -0.8848350 -0.7279159 0.31505970
6 -0.2708675 -0.9111673 -0.5143283 0.55531676
study calculated.recovery.rate
1 Arias-Gonzalez et al (2017) 4.32
2 Arias-Gonzalez et al (2017) 4.81
3 Manzello (2018) NA
4 Madin et al (2018) NA
5 Buglass et al (2015) -2.10
6 Buglass et al (2015) 0.55
# 2.) We are very curious if coral recovery rate varies across ocean basin/region. Let’s have a quick look! Make a scatterplot of recovery rate by region, if n<3 in a region, remove it. Next, make a single plot (use the ggridges package!) that shows histograms of recovery rate by region. Use patchwork to place these two graphs in the same output below.
#prelim scatter
scatterplot <- ggplot(coral_recovery_data, aes(x = region, y = calculated.recovery.rate)) +
geom_point() +
labs(y = "Recovery Rate", title = "Coral Recovery Rate by Ocean Basin/Region")
scatterplot
#filter step
coral_recovery_data2 <- coral_recovery_data %>% filter(calculated.recovery.rate >= 3)
head(coral_recovery_data2) region disturbance pre.dist.cc recovery_time hii100km
1 Caribbean Bleaching -0.22330608 -1.066311234 1.7127547
2 Caribbean Disease, Dredging -0.36202416 -1.066311234 1.7127547
3 E. Pacific COTS 0.09710549 0.001648085 0.2765956
4 E. Pacific Bleaching 0.68352937 -1.066311234 -0.5158527
5 E. Pacific Storm 0.77709608 0.802617575 -0.5158527
6 W. Pacific Bleaching, COTS 1.59961269 0.268637915 0.1799439
distance_to_shore_m dist_to_riv_m grav_NC cml_scr
1 -0.2723830 -0.8208457 -1.35285486 -0.04761585
2 -0.2723830 -0.8208457 -1.35285486 -0.04761585
3 -0.2487754 2.1269382 NA -0.21791295
4 -0.2819426 2.1864751 0.39748459 0.96539880
5 -0.2819426 2.1864751 0.39748459 0.96539880
6 0.7395376 -0.5041049 -0.02013796 -1.68293285
study calculated.recovery.rate
1 Arias-Gonzalez et al (2017) 4.320000
2 Arias-Gonzalez et al (2017) 4.810000
3 Rodgers et al (2015) 3.064000
4 Edmunds et al (2015) 6.410000
5 Edmunds et al (2015) 3.378750
6 Graham et al (2014) 3.186667#final scatter
scatterplot2 <- ggplot(coral_recovery_data2, aes(x = region, y = calculated.recovery.rate)) +
geom_point()
scatterplot2
#histogram
histogram <- ggplot(coral_recovery_data2, aes(x = calculated.recovery.rate, y = region, group = region, fill = region)) +
geom_density_ridges(stat = "binline", bins = 20, scale = 1) +
scale_y_discrete(expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
theme_ridges()
histogram
scatterplot2 / histogram
# 3.) Now, do the same for disturbance (or type of disturbance). If n<3, remove those rows/disturbances. Then remake the scatterplot and make a set of histrograms. Use patchwork to put the two graps in the same output below.
#prelim scatter
scatterplot3 <- ggplot(coral_recovery_data, aes(x = disturbance, y = calculated.recovery.rate)) +
geom_point()
scatterplot3
#filter step
# using scatterplot2 data now (which is filtered)
#final scatter
scatterplot4 <- ggplot(coral_recovery_data2, aes(x = disturbance, y = calculated.recovery.rate)) +
geom_point()
scatterplot4
#histo
histogram2 <- ggplot(coral_recovery_data2, aes(x = calculated.recovery.rate, y = disturbance, group = disturbance, fill = disturbance)) +
geom_density_ridges(stat = "binline", bins = 20, scale = 1) +
scale_y_discrete(expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
theme_ridges()
histogram2
#patchwork
scatterplot4 / histogram2
4: Our research question is– “Do local impacts effect the recovery rate of coral reefs on a global scale?” Let’s build some models to try to address this question.
a.) Does hii100km effect calculate recovery rate? Assess with an lm and a graph! Report your results.
df1 <- lm(data = coral_recovery_data2, calculated.recovery.rate ~ hii100km)
summary(df1)
Call:
lm(formula = calculated.recovery.rate ~ hii100km, data = coral_recovery_data2)
Residuals:
Min 1Q Median 3Q Max
-2.8776 -1.7345 -0.6876 0.7789 11.4392
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.3714 0.3808 14.105 <2e-16 ***
hii100km 0.3561 0.3615 0.985 0.33
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.613 on 46 degrees of freedom
Multiple R-squared: 0.02066, Adjusted R-squared: -0.0006329
F-statistic: 0.9703 on 1 and 46 DF, p-value: 0.3298plot <- ggplot(data = coral_recovery_data2, aes(x = hii100km, y = calculated.recovery.rate)) +
geom_point() +
labs(title = "HII vs Recovery Rate") +
xlab("human influence index (HII)") +
ylab("calculated recovery rate")
plot
My response HII doesn’t effect calculated recovery rate. This is clear first based on the p-value, 0.33, which is much greater than 0.05. Additionally, looking at the graph, it’s clear that there’s no significant linear relationship between HII and recovery rate.
## b.) Does cml_scr effect calculated recovery rate? Assess with an lm and a graph. Report your results.
df2 <- lm(data = coral_recovery_data2, calculated.recovery.rate ~ cml_scr)
summary(df2)
Call:
lm(formula = calculated.recovery.rate ~ cml_scr, data = coral_recovery_data2)
Residuals:
Min 1Q Median 3Q Max
-2.9512 -1.6133 -0.5723 0.8583 11.1491
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.5283 0.3907 14.148 <2e-16 ***
cml_scr 0.4721 0.3648 1.294 0.202
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.636 on 44 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.03666, Adjusted R-squared: 0.01477
F-statistic: 1.675 on 1 and 44 DF, p-value: 0.2024plot2 <- ggplot(data = coral_recovery_data2, aes(x = cml_scr, y = calculated.recovery.rate)) +
geom_point() +
labs(title = "Cumulative WCS Score vs Recovery Rate") +
xlab("Cumulative WCS Score") +
ylab("Calculated Recovery Rate")
plot2
my response Cumulative WCS score doesn’t effect calculated recovery rate. This is because the p-value of 0.202 is much greater than 0.05. Additionally, looking at the graph, it’s clear that there’s no significant linear relationship between cumulative WCS and recovery rate.
# 5.) We are really interested in more complex relationships. So, using the variables we have chosen do the following:
a.) Using a bottom-up model building structure, generate mixed effects (lmer()) models (using the lme4 pacakge) that explore the effects of the above variables on calculated recovery rate. Start with a calculated recovery rate ~ hii100km and add terms until the model starts to have singular fits. Be sure to check assumptions and for singular fit when you run each model.
We will test the following assumptions:
1. multicollinearity - essentially tests whether 2 variables are highly correlated. If they are, we struggle to determine which one is actually effecting the model. When this happens we have to make the model less complex or possibly remove a variable. We test this with check_collinearity(model_name)
2. Outliers - too many outliers will mess with our model! We can use check_outliers(model_name) for this. If we have outliers, we at least need to make a note of them.
3. A comprehensive model check for all assumptions (normality of residuals, homogeneity of variance, outliers, collinearity). We use check_model(model_name) from the performance package for this.
IF a model has a singular fit or fails the collinearity assumption, we cannot use it. You either have too complex a model or variables that are collinear. You can always regress variables or use a correlation plot to see if you have collinear variables.
ALL lmer models MUST have a random effect. We will use study as our random effect. You add this to the model by doing adding + (1|study) in the model.
While doing model testing we must use the option REML=FALSE in the model. REML is reduced maximum likelihood. We will change REML to TRUE once we select a final model.
Fixed vs. Random effects: Fixed don’t vary across groups or the effect is ‘fixed’. Random = account for vars that differ between individuals. Also something that is a covariate that you aren’t expecting to be part of the design / didn’t plan for.
#model 1 should be
lmer1<-lmer(calculated.recovery.rate~hii100km + (1|study), data = coral_recovery_data2, REML=FALSE)
lmer1Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km + (1 | study)
Data: coral_recovery_data2
AIC BIC logLik deviance df.resid
232.9300 240.4148 -112.4650 224.9300 44
Random effects:
Groups Name Std.Dev.
study (Intercept) 1.444
Residual 2.131
Number of obs: 48, groups: study, 31
Fixed Effects:
(Intercept) hii100km
5.3074 0.3276 check_collinearity(lmer1)Not enough model terms in the conditional part of the model to check for
multicollinearity.NULLcheck_outliers(lmer1)Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.2 outliers detected: cases 11, 43.
- Based on the following method and threshold: cook (0.704).
- For variable: (Whole model).check_model(lmer1)Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.
#model 2 should be y*x
lmer2<-lmer(calculated.recovery.rate~hii100km*cml_scr + (1|study), data = coral_recovery_data2, REML=FALSE)
df2 <- lm(data = coral_recovery_data2, calculated.recovery.rate~hii100km*cml_scr)
lmer2Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km * cml_scr + (1 | study)
Data: coral_recovery_data2
AIC BIC logLik deviance df.resid
227.6568 238.6287 -107.8284 215.6568 40
Random effects:
Groups Name Std.Dev.
study (Intercept) 1.428
Residual 2.141
Number of obs: 46, groups: study, 30
Fixed Effects:
(Intercept) hii100km cml_scr hii100km:cml_scr
5.3599 -0.2287 0.5328 0.2666 check_collinearity(lmer2)# Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
hii100km 2.45 [1.74, 3.84] 1.56 0.41 [0.26, 0.57]
cml_scr 1.46 [1.16, 2.30] 1.21 0.69 [0.43, 0.86]
hii100km:cml_scr 2.03 [1.49, 3.16] 1.42 0.49 [0.32, 0.67]check_outliers(lmer2) Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.2 outliers detected: cases 11, 41.
- Based on the following method and threshold: cook (0.801).
- For variable: (Whole model).check_model(df2) Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.
Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.
- region (what ocean basin is the coral from?)
- disturbance (type of disturbance that occurred)
- pre.dist.cc (the coral cover at the site before the disturbance in %)
- recovery_time (time in years over which recovery rate is measured)
- hii100km (the human influence index in a 100km radius of the site)
- distance_to_shore_m (the distance to the nearest major shoreline in meters)
- dist_to_riv_m (the distance to the nearest major river mount in meters)
- grav_NC (market gravity - proximity to fish market. This is actually the number of people within a 5km range of the reef)
- cml_scr (cumulative WCS score– a cumulative measure of human influence - similar to HII, but updated to be reef specific)
- study (the published paper that these data came from)
#model 3 should be y~x*z
lmer3 <-lmer(calculated.recovery.rate~hii100km*cml_scr*recovery_time + disturbance + (1|study), data = coral_recovery_data2, REML=FALSE)
lmer3Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km * cml_scr * recovery_time +
disturbance + (1 | study)
Data: coral_recovery_data2
AIC BIC logLik deviance df.resid
229.5284 256.9580 -99.7642 199.5284 31
Random effects:
Groups Name Std.Dev.
study (Intercept) 0.8295
Residual 1.9588
Number of obs: 46, groups: study, 30
Fixed Effects:
(Intercept) hii100km
5.14948 -0.63141
cml_scr recovery_time
0.79744 -0.68930
disturbanceBleaching, COTS disturbanceBleaching, Storm
-0.46011 0.28763
disturbanceCOTS disturbanceDisease, Dredging
1.55860 0.12567
disturbanceStorm hii100km:cml_scr
-0.30333 0.54934
hii100km:recovery_time cml_scr:recovery_time
-0.05567 -0.59211
hii100km:cml_scr:recovery_time
-0.48539 check_collinearity(lmer3)# Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance
hii100km 3.68 [2.69, 5.26] 1.92 0.27
cml_scr 1.87 [1.47, 2.61] 1.37 0.54
recovery_time 2.03 [1.58, 2.84] 1.43 0.49
disturbance 3.71 [2.71, 5.29] 1.93 0.27
hii100km:cml_scr 2.86 [2.14, 4.05] 1.69 0.35
hii100km:recovery_time 3.69 [2.69, 5.27] 1.92 0.27
cml_scr:recovery_time 1.70 [1.36, 2.36] 1.30 0.59
hii100km:cml_scr:recovery_time 3.64 [2.66, 5.20] 1.91 0.27
Tolerance 95% CI
[0.19, 0.37]
[0.38, 0.68]
[0.35, 0.63]
[0.19, 0.37]
[0.25, 0.47]
[0.19, 0.37]
[0.42, 0.74]
[0.19, 0.38]check_outliers(lmer3) Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.OK: No outliers detected.
- Based on the following method and threshold: cook (0.9).
- For variable: (Whole model)check_model(lmer3) Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.
Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.
#model 4 should be y
lmer4<-lmer(calculated.recovery.rate ~ hii100km*cml_scr*region + (1|study), data = coral_recovery_data2, REML=FALSE)boundary (singular) fit: see help('isSingular')lmer4Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km * cml_scr * region + (1 |
study)
Data: coral_recovery_data2
AIC BIC logLik deviance df.resid
240.8800 273.7956 -102.4400 204.8800 28
Random effects:
Groups Name Std.Dev.
study (Intercept) 0.000
Residual 2.243
Number of obs: 46, groups: study, 30
Fixed Effects:
(Intercept) hii100km
6.5179 -1.0379
cml_scr regionE. Pacific
3.2473 -5.7467
regionIndian Ocean regionW. Pacific
-2.3280 0.4842
hii100km:cml_scr hii100km:regionE. Pacific
-0.7166 8.7754
hii100km:regionIndian Ocean hii100km:regionW. Pacific
-0.7291 1.4165
cml_scr:regionE. Pacific cml_scr:regionIndian Ocean
-0.6832 -2.2819
cml_scr:regionW. Pacific hii100km:cml_scr:regionE. Pacific
-1.6897 -11.0861
hii100km:cml_scr:regionIndian Ocean hii100km:cml_scr:regionW. Pacific
1.8814 -0.0362
optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings check_collinearity(lmer4)Model has interaction terms. VIFs might be inflated.
You may check multicollinearity among predictors of a model without
interaction terms.# Check for Multicollinearity
High Correlation
Term VIF VIF 95% CI Increased SE Tolerance
hii100km 53.93 [ 38.06, 76.60] 7.34 0.02
cml_scr 153.96 [ 108.33, 218.99] 12.41 6.50e-03
region 1047.94 [ 736.35, 1491.58] 32.37 9.54e-04
hii100km:cml_scr 35.68 [ 25.24, 50.63] 5.97 0.03
hii100km:region 22567.44 [15853.63, 32124.63] 150.22 4.43e-05
cml_scr:region 19450.80 [13664.21, 27688.09] 139.47 5.14e-05
hii100km:cml_scr:region 10291.33 [ 7229.76, 14649.56] 101.45 9.72e-05
Tolerance 95% CI
[0.01, 0.03]
[0.00, 0.01]
[0.00, 0.00]
[0.02, 0.04]
[0.00, 0.00]
[0.00, 0.00]
[0.00, 0.00]check_outliers(lmer4)Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.3 outliers detected: cases 7, 25, 26.
- Based on the following method and threshold: cook (0.801).
- For variable: (Whole model).check_model(lmer4)Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.Some of the variables were in matrix-format - probably you used
`scale()` on your data?
If so, and you get an error, please try `datawizard::standardize()` to
standardize your data.
b.) Check the performance of ALL of the models that did not result in an error (a singular fit) to see which one is the best fit for the data. We can use the compare_performance() function from the performance package for this. The syntax is compare_performance(lm1,lm2,lm3, rank=TRUE). This rank=True will put the ‘best’ model at the top.
Compare performance will print some combination of the following metrics:
1. R2 and R2 adjusted - this is a great metric that we already understand (goodness of fit!)
2. RMSE - root mean squared error (smaller is generally better) 3. AIC, BIC are both criterion for model selection-> smaller is better. If the column says AIC weights, BIC weights, or AICc weights, then bigger is better! 4. ICC is Intraclass correlation coefficient - “the proportion of the variance explained by the grouping structure in the population” (Hox 2010) 5. BF is bayes factor (not sure it is relevant) 6. Finally- we often get a Performance-Score, which is an estimated goodness of fit value from the package. This is NOT ALWAYS the best way to choose a model. BUT, it often tells us at least which model is the best fit (statistically). Notably, stats are a tool– if the model doesn’t address our hypothesis, we will often not use the “best” model.
#compare model performance
compare_performance(lmer1,lmer2,lmer3, rank=TRUE)When comparing models, please note that probably not all models were fit
from same data.# Comparison of Model Performance Indices
Name | Model | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | AIC weights | AICc weights | BIC weights | Performance-Score
--------------------------------------------------------------------------------------------------------------------------------
lmer3 | lmerMod | 0.445 | 0.346 | 0.152 | 1.821 | 1.959 | 0.268 | 3.41e-04 | 7.43e-05 | 54.32%
lmer2 | lmerMod | 0.345 | 0.054 | 0.308 | 1.853 | 2.141 | 0.683 | 0.883 | 0.709 | 52.84%
lmer1 | lmerMod | 0.327 | 0.018 | 0.315 | 1.839 | 2.131 | 0.049 | 0.117 | 0.290 | 25.57%my best fitting model is lmer3, which is tests the relationship, calculated.recovery.rate~hii100km * cml_scr * recovery_time + (1|study), is the best fitting model. It has the highest performance score and the values for both conditional R^2 and margincal R^2 are the largest.
Now, add REML=TRUE to your best model and compare it to the same model with REML=FALSE. Notice any difference?
summary(lmer3)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km * cml_scr * recovery_time +
disturbance + (1 | study)
Data: coral_recovery_data2
AIC BIC logLik deviance df.resid
229.5 257.0 -99.8 199.5 31
Scaled residuals:
Min 1Q Median 3Q Max
-1.4119 -0.6100 -0.0149 0.1810 3.9131
Random effects:
Groups Name Variance Std.Dev.
study (Intercept) 0.6881 0.8295
Residual 3.8371 1.9588
Number of obs: 46, groups: study, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.14948 0.63106 8.160
hii100km -0.63141 0.59246 -1.066
cml_scr 0.79744 0.42041 1.897
recovery_time -0.68930 0.42717 -1.614
disturbanceBleaching, COTS -0.46011 2.38558 -0.193
disturbanceBleaching, Storm 0.28763 2.21924 0.130
disturbanceCOTS 1.55860 1.21454 1.283
disturbanceDisease, Dredging 0.12567 2.47127 0.051
disturbanceStorm -0.30333 0.98107 -0.309
hii100km:cml_scr 0.54934 0.43915 1.251
hii100km:recovery_time -0.05567 0.60200 -0.092
cml_scr:recovery_time -0.59211 0.37195 -1.592
hii100km:cml_scr:recovery_time -0.48539 0.47608 -1.020
Correlation matrix not shown by default, as p = 13 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need itlmer3_reml_false <-lmer(calculated.recovery.rate~hii100km*cml_scr*recovery_time + (1|study), data = coral_recovery_data2, REML=TRUE)
summary(lmer3_reml_false)Linear mixed model fit by REML ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km * cml_scr * recovery_time +
(1 | study)
Data: coral_recovery_data2
REML criterion at convergence: 203.1
Scaled residuals:
Min 1Q Median 3Q Max
-1.0454 -0.5711 -0.0953 0.1557 3.2681
Random effects:
Groups Name Variance Std.Dev.
study (Intercept) 1.797 1.341
Residual 4.257 2.063
Number of obs: 46, groups: study, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.3713 0.4798 11.195
hii100km -0.3582 0.6166 -0.581
cml_scr 0.7593 0.4584 1.657
recovery_time -0.4990 0.4102 -1.217
hii100km:cml_scr 0.3497 0.4670 0.749
hii100km:recovery_time 0.2755 0.6536 0.422
cml_scr:recovery_time -0.6741 0.4146 -1.626
hii100km:cml_scr:recovery_time -0.6643 0.5316 -1.250
Correlation of Fixed Effects:
(Intr) h100km cml_sc rcvry_ h100km:c_ h100km:r_ cml_:_
hii100km 0.095
cml_scr 0.177 -0.386
recovery_tm -0.140 0.093 0.054
h100km:cml_ -0.434 -0.686 -0.008 -0.044
h100km:rcv_ 0.125 0.282 0.118 0.058 -0.300
cml_scr:rc_ 0.090 0.187 -0.214 0.327 -0.157 -0.203
h100km:c_:_ -0.083 -0.296 -0.095 -0.381 0.281 -0.730 -0.176The main difference is that the disturbance variable is not measured when we set REML=TRUE. Additionally, the standard deviations of the random effects are larger for REML=TRUE and the interactions fixed effects values changed as well.
c.) Make a 95% CI coefficients table for your model. Tell me which of your model parameters/variables have positive/negative effects. Examples of this kind of plot can be found here
coefs<-tidy(lmer3, quick=FALSE)
ci<-data.table(confint(lmer3), keep.rownames='term')Computing profile confidence intervals ...cidf<-cbind(coefs,ci)
cidf<-cidf[,-6]
cidf<- cidf %>%
rename("lower"="2.5 %",
"upper"="97.5 %")
cidf$term=as.factor(cidf$term)
cidf effect group term estimate std.error
1 fixed <NA> (Intercept) 5.14948030 0.6310609
2 fixed <NA> hii100km -0.63140817 0.5924618
3 fixed <NA> cml_scr 0.79744347 0.4204087
4 fixed <NA> recovery_time -0.68929976 0.4271701
5 fixed <NA> disturbanceBleaching, COTS -0.46011258 2.3855777
6 fixed <NA> disturbanceBleaching, Storm 0.28763360 2.2192414
7 fixed <NA> disturbanceCOTS 1.55860220 1.2145423
8 fixed <NA> disturbanceDisease, Dredging 0.12567487 2.4712733
9 fixed <NA> disturbanceStorm -0.30333158 0.9810737
10 fixed <NA> hii100km:cml_scr 0.54933970 0.4391475
11 fixed <NA> hii100km:recovery_time -0.05566584 0.6019957
12 fixed <NA> cml_scr:recovery_time -0.59211462 0.3719503
13 fixed <NA> hii100km:cml_scr:recovery_time -0.48538647 0.4760757
14 ran_pars study sd__(Intercept) 0.82951943 NA
15 ran_pars Residual sd__Observation 1.95884794 NA
term.1 lower upper
1 .sig01 0.00000000 1.8681195
2 .sigma 1.44430170 2.6204623
3 (Intercept) 3.85738691 6.4130883
4 hii100km -1.84427610 0.5545549
5 cml_scr -0.04850728 1.6937212
6 recovery_time -1.54538908 0.1658695
7 disturbanceBleaching, COTS -5.24040776 4.3738729
8 disturbanceBleaching, Storm -4.15481984 4.7404870
9 disturbanceCOTS -0.89601247 4.0455863
10 disturbanceDisease, Dredging -5.14215385 5.0722518
11 disturbanceStorm -2.26709142 1.6902678
12 hii100km:cml_scr -0.34147067 1.4294338
13 hii100km:recovery_time -1.26067291 1.2003137
14 cml_scr:recovery_time -1.33754796 0.1529004
15 hii100km:cml_scr:recovery_time -1.45997778 0.4675923ggplot(data=cidf, aes(x=estimate, y=term.1))+
geom_vline(xintercept = 0, linetype=2)+
geom_errorbarh(aes(xmax=lower, xmin=upper),height=0.2)+
theme_classic()
The cml_scr variable has a positive effect. All others don’t have significant positive or negative effects because their 95% confidence intervals all lie within the dotted line. You could also say that cml_scr:recovery_time most likely has a negative effect because while it does lie within the dotted line, it is almost completely to the right of it.
d.) Make plots to assess the effect of your explanatory variables. You will need to extract the model fits for each var and plot!
This requires ggpredict() from the ggeffects package! ggpredict produces 95% CI from the model! This is really useful for visual assessment of results.
Some example code:
#run the linear model
#finalmodel2<-lmer(y~x*z*a)
# Extract the prediction data frame
#pred.hii <- ggpredict(finalmodel2, terms = c("hii100km2")) # this gives overall predictions for the model
#pred.hii
#
# # Plot the predictions
# ggpredhii<-ggplot(pred.hii) +
# geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high),
# fill = "lightgrey", alpha = 0.3) + # error band
# geom_point(data = recov, # adding the raw data (scaled values)
# aes(x = hii100km2, y = calculated.recovery.rate),alpha=0.5) +
# geom_line(aes(x = x, y = predicted),size=1, color='black') + # slope
# theme_bw()+
# labs(x = "Scaled Human Influence Index", y = "Recovery Rate")
## AND A CATEGORICAL VAR
###REGION
# pred.region <- ggpredict(finalmodel2, terms = c("region")) # this gives overall predictions for the model
# summary(pred.region)
# pred.region
# pred.region$x<-factor(pred.region$x, levels= c("Caribbean","Indian Ocean", "W. Pacific", "E. Pacific"))
#
# ggpredreg<-ggplot(pred.region) +
# geom_point(data=recov, aes(x=region, y= calculated.recovery.rate), color='grey', alpha=0.3)+
# geom_point(aes(x = x, y = predicted, colour = x),size=3) +
# geom_errorbar(aes(x=x, ymin=conf.low, ymax=conf.high, color=x),width=0.5)+
# #geom_errorbar(aes(x=x, ymin=predicted-std.error, ymax=predicted+std.error),width=0.3)+
# theme_bw()+
# labs(x = "Region", y = "Recovery Rate")
YOU TRY - Do this for a least 1 numerical and 1 categorical explanatory variable. Interpret your results based on your model summary table, coefficient plot, and these graphs. What do you see?
#run the linear model
#finalmodel2<-lmer(y~x*z*a)
lmer3Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: calculated.recovery.rate ~ hii100km * cml_scr * recovery_time +
disturbance + (1 | study)
Data: coral_recovery_data2
AIC BIC logLik deviance df.resid
229.5284 256.9580 -99.7642 199.5284 31
Random effects:
Groups Name Std.Dev.
study (Intercept) 0.8295
Residual 1.9588
Number of obs: 46, groups: study, 30
Fixed Effects:
(Intercept) hii100km
5.14948 -0.63141
cml_scr recovery_time
0.79744 -0.68930
disturbanceBleaching, COTS disturbanceBleaching, Storm
-0.46011 0.28763
disturbanceCOTS disturbanceDisease, Dredging
1.55860 0.12567
disturbanceStorm hii100km:cml_scr
-0.30333 0.54934
hii100km:recovery_time cml_scr:recovery_time
-0.05567 -0.59211
hii100km:cml_scr:recovery_time
-0.48539 # Extract the prediction data frame
pred.hii <- ggpredict(lmer3, terms = c("hii100km")) # this gives overall predictions for the model
pred.hii# Predicted values of calculated.recovery.rate
hii100km | Predicted | 95% CI
----------------------------------
-0.80 | 5.49 | 3.95, 7.03
-0.40 | 5.21 | 3.87, 6.55
0.20 | 4.80 | 3.36, 6.23
0.60 | 4.52 | 2.78, 6.26
1.00 | 4.24 | 2.11, 6.38
1.60 | 3.83 | 0.99, 6.66
2.00 | 3.55 | 0.22, 6.89
3.00 | 2.86 | -1.77, 7.49
Adjusted for:
* cml_scr = -0.11
* recovery_time = 0.20
* disturbance = Bleaching
* study = 0 (population-level)
Not all rows are shown in the output. Use `print(..., n = Inf)` to show
all rows.#
# # Plot the predictions
ggpredhii<-ggplot(pred.hii) +
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), fill = "lightgrey", alpha = 0.3) + # error band
geom_point(data = coral_recovery_data2, # adding the raw data (scaled values)
aes(x = hii100km, y = calculated.recovery.rate),alpha=0.5) +
geom_line(aes(x = x, y = predicted),size=1, color='black') + # slope
theme_bw()+
labs(x = "Scaled Human Influence Index", y = "Recovery Rate")Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.ggpredhii
## AND A CATEGORICAL VAR
pred.disturbance <- ggpredict(lmer3, terms = c("disturbance")) # this gives overall predictions for the model
summary(pred.disturbance) x predicted std.error conf.low
Bleaching :1 Min. :4.347 Min. :0.7008 Min. :-0.2001
Bleaching, COTS :1 1st Qu.:4.580 1st Qu.:0.7981 1st Qu.: 0.1483
Bleaching, Storm :1 Median :4.870 Median :1.5454 Median : 1.7735
COTS :1 Mean :5.009 Mean :1.5349 Mean : 1.8782
Disease, Dredging:1 3rd Qu.:5.055 3rd Qu.:2.2223 3rd Qu.: 3.2683
Storm :1 Max. :6.366 Max. :2.4211 Max. : 4.5494
conf.high group
Min. :6.069 1:6
1st Qu.:6.723
Median :8.539
Mean :8.139
3rd Qu.:9.410
Max. :9.871 pred.disturbance# Predicted values of calculated.recovery.rate
disturbance | Predicted | 95% CI
-------------------------------------------
Bleaching | 4.81 | 3.38, 6.24
Bleaching, COTS | 4.35 | -0.20, 8.89
Bleaching, Storm | 5.10 | 0.61, 9.58
COTS | 6.37 | 4.55, 8.18
Disease, Dredging | 4.93 | 0.00, 9.87
Storm | 4.50 | 2.94, 6.07
Adjusted for:
* hii100km = 0.18
* cml_scr = -0.11
* recovery_time = 0.20
* study = 0 (population-level)pred.disturbance$x<-factor(pred.disturbance$x, levels= c("Storm","Disease, Dredging", "COTS", "Bleaching, Storm", "Bleaching, COTS"))
#
ggpredreg<-ggplot(pred.disturbance) +
geom_point(data=coral_recovery_data2, aes(x=disturbance, y= calculated.recovery.rate, na.rm = TRUE), color='grey', alpha=0.3)+
geom_point(aes(x = x, y = predicted, colour = x),size=3) +
geom_errorbar(aes(x=x, ymin=conf.low, ymax=conf.high, color=x),width=0.5)+
geom_errorbar(aes(x=x, ymin=predicted-std.error, ymax=predicted+std.error),width=0.3)+
theme_bw()+
labs(x = "Disturbance", y = "Recovery Rate") +
scale_x_discrete(labels = label_wrap(10)) Warning in geom_point(data = coral_recovery_data2, aes(x = disturbance, :
Ignoring unknown aesthetics: na.rmggpredreg
My interpretation
First looking at the graphs in this section, There seems to be a slightly negative impact on recovery rate as HII increases but it’s not clear from the graph alone if it has a significant effect on recovery rate. Then looking at the model summary, specific to HII, an increase of 1 percent recovery rate would result in an average -0.63141 decrease in hii. Finally, looking at the confidence interval plot, it doesn’t appear that HII has a significant impact on recovery rate. Thus, HII doesn’t have a significant influence on recovery rate. Now looking at the graphs in this section for disturbance, the average recovery doesn’t seem to vary much by type of disturbance. However, there is a greater variability in rate for disease and bleaching related disturbances than for COTS and storm disturbances. For the model summary, the disturbances have varied effects of increasing or decreasing average recovery rate, with COTS having the greatest effect. Then finally looking at the confidence interval plot, none of the disturbances significantly varied the average recovery rate, meaning that on it’s own disturbance doesn’t have a significant effect on recovery rate.
# Render the lab and turn it in :)