ANLY 505 - Problem Set #2
Simulation of Linear Mixed Models
Linto Lalu
2019-10-14
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 = 1Simulate 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):
- 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)## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 3.5.3fit2 = lmer(resp2 ~ elevation + slope + (1|stand))
fit2## Linear mixed model fit by REML ['lmerMod']
## Formula: resp2 ~ 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.09511summary(fit2)## Linear mixed model fit by REML ['lmerMod']
## Formula: resp2 ~ 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 t value
## (Intercept) -21.314628 6.602050 -3.228
## elevation 0.020600 0.004916 4.190
## slope 0.095105 0.016441 5.785
##
## Correlation of Fixed Effects:
## (Intr) elevtn
## elevation -0.991
## slope 0.049 -0.148The b0 estimate= -21.314628. Widely away far from -1. The estimate of b1 is 0.0206 that is not close to .005 but estimate of b2 is 0.095 is close to 0.1. The estimate of standard deviation for stand is 1.099
- Create a function for your model and run 1000 simulations of that model.
set.seed(16)
fit2model_fun = function(nstand = 5, nplot = 4, b0 = -1, b1 = 0.005, b2 = .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
dframe = data.frame(resp2, elevation, slope, stand)
fit2 = lmer(resp2 ~ elevation + slope + (1|stand), data = dframe)
}
set.seed(16)
sim_fit2model = replicate(1000 , expr = fit2model_fun())## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit- Extract the stand and residual variances from this simulation run. Print the first 6 rows of the data.
library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.5.3## -- Attaching packages ----------------------------------------------- tidyverse 1.2.1 --## v ggplot2 3.2.0 v purrr 0.3.2
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.3.1
## v readr 1.3.1 v forcats 0.3.0## Warning: package 'ggplot2' was built under R version 3.5.3## Warning: package 'tibble' was built under R version 3.5.3## Warning: package 'tidyr' was built under R version 3.5.3## Warning: package 'purrr' was built under R version 3.5.3## Warning: package 'dplyr' was built under R version 3.5.3## -- Conflicts -------------------------------------------------- tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()library(broom)## Warning: package 'broom' was built under R version 3.5.3library(purrr)
sv = sim_fit2model %>%
map_dfr(tidy, effects = "ran_pars",scales = "vcov")
head(sv)## # A tibble: 6 x 3
## term group estimate
## <chr> <chr> <dbl>
## 1 var_(Intercept).stand stand 1.21
## 2 var_Observation.Residual Residual 1.36
## 3 var_(Intercept).stand stand 5.56
## 4 var_Observation.Residual Residual 0.951
## 5 var_(Intercept).stand stand 2.61
## 6 var_Observation.Residual Residual 1.11- 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?
stand_sims = c(5, 20, 100) %>%
set_names(c("sample size 5","sample size 20","sample size 100")) %>%
map(~replicate(1000, fit2model_fun(nstand = .x)))## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fit
## singular fitlibrary(broom)
stand_vars = stand_sims %>%
modify_depth(2, ~tidy(.x, effects = "ran_pars", scales = "vcov")) %>%
map_dfr(bind_rows, .id = "stand_num") %>%
filter(group == "stand")
stand_vars = mutate(stand_vars, stand_num = fct_inorder(stand_num))
add_prefix = function(string) {
paste("# of stands:" , string, sep = " ")
}
groupmed = stand_vars %>%
group_by(stand_num) %>%
summarise(mvar = median(estimate))
library(ggplot2)
ggplot(stand_vars, aes(x = estimate))+
geom_density(fill = "orange", alpha = 0.25) +
facet_wrap(~stand_num, labeller = as_labeller(add_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() +
scale_linetype_manual(name = "", values = c(2,1)) +
theme(legend.position = "bottom", legend.key.width = unit(.1, "cm")) +
labs(x = "Estimated Variance", y = "Density")
- Plot the coefficients of the estimates of elevation and slope. Hint: the x-axis should have 1000 values. Discuss the graphs.
coefplot = sim_fit2model %>%
map_dfr(tidy, effects = "fixed") %>%
filter(term %in% c("elevation","slope")) %>%
group_by(term)
coefplot## # A tibble: 2,000 x 4
## # Groups: term [2]
## term estimate std.error statistic
## <chr> <dbl> <dbl> <dbl>
## 1 elevation 0.0206 0.00492 4.19
## 2 slope 0.0951 0.0164 5.78
## 3 elevation -0.00546 0.00927 -0.589
## 4 slope 0.0868 0.0140 6.21
## 5 elevation 0.0152 0.00426 3.57
## 6 slope 0.121 0.0157 7.73
## 7 elevation 0.0156 0.0117 1.33
## 8 slope 0.124 0.0166 7.48
## 9 elevation 0.0130 0.00336 3.87
## 10 slope 0.0937 0.0100 9.32
## # ... with 1,990 more rowslibrary(dplyr)
mean = coefplot %>%
summarise_at("estimate", funs( mean))## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.mean## # A tibble: 2 x 2
## term estimate
## <chr> <dbl>
## 1 elevation 0.00492
## 2 slope 0.0997plot = ggplot(coefplot, aes(x=rep(1:1000, each=2),y=estimate))+
geom_point(size=0.5)+
facet_wrap(~term) +
geom_hline(data=mean,aes(yintercept=mean$estimate,color = term), size = 0.9) +
labs(x="Replicated numbers", y="Estimate of coefficients")
plot