The statistical model:

\(y_t = \beta_0 + \beta_1 * (Elevation_s)_t + \beta_2 * Slope_t + (b_s)_t + \epsilon_t\)

Where:

• \(\beta_0\) is the mean response when both Elevation and Slope are 0
• \(\beta_1\) is the change in mean response for a 1-unit change in elevation. Elevation is measured at the stand level, so all plots in a stand share a single value in elevation.
• \(\beta_2\) is the change in mean response for a 1-unit change in slope. Slope is measured at the plot level, so every plot potentially has a unique value of slope.

Let’s define the parameters:

• the intercept \((\beta_0)\) will be -1
• the coefficient for elevation \((\beta_1)\) will be set to 0.005
• the coefficient for slope \((\beta_2)\) will be set to 0.1

nstand = 5
nplot = 4
b0 = -1
b1 = .005
b2 = .1
sds = 2
sd = 1

Simulate other variables:

set.seed(16)
stand = rep(LETTERS[1:nstand], each = nplot)
standeff = rep( rnorm(nstand, 0, sds), each = nplot)
ploteff = rnorm(nstand*nplot, 0, sd)

Simulate elevation and slope:

elevation = rep( runif(nstand, 1000, 1500), each = nplot)
slope = runif(nstand*nplot, 2, 75)

Simulate response variable:

resp2 = b0 + b1*elevation + b2*slope + standeff + ploteff 

Your tasks (complete each task in its’ own code chunk, make sure to use echo=TRUE so I can see your code):

1. Fit a linear mixed model with the response variable as a function of elevation and slope with stand as a random effect. Are the estimated parameters similar to the intial parameters as we defined them?

library(lme4)

library(Matrix)

m1fit = lmer(resp2 ~ 1 + elevation + slope + (1|stand))

m1fit

Linear mixed model fit by REML [‘lmerModLmerTest’] Formula: resp2 ~ 1 + elevation + slope + (1 | stand) REML criterion at convergence: 81.9874 Random effects: Groups Name Std.Dev. stand (Intercept) 1.099
Residual 1.165
Number of obs: 20, groups: stand, 5 Fixed Effects: (Intercept) elevation slope
-21.31463 0.02060 0.09511##

summary(m1fit) Linear mixed model fit by REML. t-tests use Satterthwaite’s method [ lmerModLmerTest] Formula: resp2 ~ 1 + elevation + slope + (1 | stand)

REML criterion at convergence: 82

Scaled residuals: Min 1Q Median 3Q Max -1.65583 -0.62467 -0.01693 0.53669 1.41736

Random effects: Groups Name Variance Std.Dev. stand (Intercept) 1.208 1.099
Residual 1.358 1.165
Number of obs: 20, groups: stand, 5

Fixed effects: Estimate Std. Error df t value Pr(>|t|)
(Intercept) -21.314628 6.602053 3.001313 -3.228 0.0482 *
elevation 0.020600 0.004916 3.113482 4.190 0.0230 *
slope 0.095105 0.016441 15.868032 5.785 2.88e-05 *** — Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects: (Intr) elevtn elevation -0.991
slope 0.049 -0.148

2. Create a function for your model and run 1000 simulations of that model.

mixed_function = function(nstand = 5, nplot = 4, b0 = -1, b1 = 0.005, b2 = 0.1, sds = 2, sd = 1) { + stand = rep(LETTERS[1:nstand], each = nplot) + standeff = rep(rnorm(nstand, 0, sds), each = nplot) + ploteff = rnorm(nstand * nplot, 0, sd) + elevation = rep(runif(nstand, 1000, 1500), each = nplot) + slope = runif(nstand * nplot, 2, 75) + resp2 = b0 + b1 * elevation + b2 * slope + standeff + ploteff + dat = data.frame(resp2, elevation, slope, stand) + lmer(resp2 ~ 1 + elevation + slope + (1|stand), data = dat) + } mixed_function()

Linear mixed model fit by REML [‘lmerModLmerTest’] Formula: resp2 ~ 1 + elevation + slope + (1 | stand) Data: dat REML criterion at convergence: 84.9965 Random effects: Groups Name Std.Dev. stand (Intercept) 1.917
Residual 1.251
Number of obs: 20, groups: stand, 5 Fixed Effects: (Intercept) elevation slope
-20.64450 0.02083 0.09546
convergence code 0; 1 optimizer warnings; 0 lme4 warnings Warning message:

repeat_function = replicate( n = 1000, expr = mixed_function())

boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular boundary (singular) fit: see ?isSingular

Warning messages: 1: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00406233 (tol = 0.002, component 1)

2: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.0032079 (tol = 0.002, component 1)

3: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00561834 (tol = 0.002, component 1)

4: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00438273 (tol = 0.002, component 1)

5: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00230691 (tol = 0.002, component 1)

6: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00280278 (tol = 0.002, component 1)

7: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00651067 (tol = 0.002, component 1)

8: In checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00555055 (tol = 0.002, component 1)

3. Extract the stand and residual variances from this simulation run. Print the first 6 rows of the data.

library(broom) variances <- simsmix %>% map_dfr(tidy, effects = “ran_pars”, scales = “vcov”) variances %>% print(n = 6)

A tibble: 2,000 x 4 effect group term estimate 1 ran_pars stand var__(Intercept) 0.966 2 ran_pars Residual varObservation 1.20
3 ran_pars stand var(Intercept) 2.91
4 ran_pars Residual varObservation 1.49
5 ran_pars stand var(Intercept) 0.0990 6 ran_pars Residual var__Observation 1.58

… with 1,994 more rows

4. Choose three different sample sizes (your choice) and run 1000 model simulations with each sample size. Create 3 visualizations that compare distributions of the variances for each of the 3 sample sizes. Make sure that the axes are labelled correctly. What do these graphs say about the relationship between sample size and variance?

library(ggplot2) library(dplyr)

simsample_size = c(10, 25, 35) %>% + set_names(c(“5”, “15”, “25”)) %>% + map(~replicate(1000, mixed_function(nstand = .x)))

boundary (singular) fit: see ?isSingular Warning in checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00446168 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, “derivs”), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.00617958 (tol = 0.002, component 1) 5. Plot the coefficients of the estimates of elevation and slope. Hint: the x-axis should have 1000 values. Discuss the graphs.

vars = simsample_size %>% + modify_depth(2, ~tidy(.x, effects = “ran_pars”, scales = “vcov”)) %>% + map_dfr(bind_rows, .id = “id”) %>% + filter(group == “stand”) prefix = function(string) { + paste(“Number of stands:”, string, sep = " ") + }

groupmed = vars %>% + group_by(id) %>% + summarise(mvar = median(estimate))

ggplot(vars, aes(x = estimate)) + + geom_density(fill = “green”, alpha = “0.25”) + + facet_wrap(~id, labeller = as_labeller(prefix)) + + geom_vline(aes(xintercept = 4, linetype = “True Variance”), size=0.5) + + geom_vline(data = groupmed, aes(xintercept=mvar,linetype=“Median Variance”), size = 0.5) + + theme_bw() + + labs(x=“Estimated Variance”)