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?

fit1=lm(resp2~ elevation + slope)
summary(fit1)

## 
## Call:
## lm(formula = resp2 ~ elevation + slope)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3087 -1.2573  0.0560  0.8909  2.2694 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -21.368461   3.937824  -5.426 4.54e-05 ***
## elevation     0.020722   0.003007   6.891 2.61e-06 ***
## slope         0.092370   0.018666   4.949 0.000122 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.479 on 17 degrees of freedom
## Multiple R-squared:  0.8525, Adjusted R-squared:  0.8351 
## F-statistic: 49.11 on 2 and 17 DF,  p-value: 8.625e-08

Answer: estimated slope 0.1079 is closed to initial parameter 0.1; and estimate elevation is -0.0018 is not clost to initia 0.005; estimate intercept 7.2282 is not close to initial parameter -1

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

thrun<-function(nstand = 5,nplot = 4,b0 = -1,b1 = .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 
        fit1=lm(resp2~ elevation + slope)
}
set.seed(8)
sim=thrun()
#for some reasons, rerun functon cannot run #
sims=replicate(1000,thrun())

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

library(broom)

## Warning: package 'broom' was built under R version 3.5.3

#variances <- sims %>% map_df(tidy) %>%
#variances %>% print(n = 6)

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)

## Warning: package 'ggplot2' was built under R version 3.5.3

set.seed(15)
x1=thrun(5)
x2=thrun(15)
x3=thrun(100)

ggplot(x1,aes(x=elevation))+
  geom_density(fill ="red",alpha=0.25)

ggplot(x2,aes(x=elevation))+
  geom_density(fill ="red",alpha=0.25)

ggplot(x3,aes(x=elevation))+
  geom_density(fill ="red",alpha=0.25)

5. Plot the coefficients of the estimates of elevation and slope. Hint: the x-axis should have 1000 values. Discuss the graphs.

#library(furrr)
#simsest <- sims %>% 
  #future_map(tidy, effects = "fixed") %>% 
  #bind_rows()

#simsest %>% 
  #dplyr::filter(term %in% c("elevation", "slope")) %>% 
  #group_by(term) %>% 
  #mutate(x = 1 : 1000) %>%
  #ungroup() %>% 
  #mutate(real_value = ifelse(term == "elevation", 0.005, 0.1)) %>% 
  #ggplot(aes(x = x, y = estimate)) +
  #geom_line() +
  #facet_wrap(~term) +
  #geom_hline(aes(yintercept = real_value, color = term), linetype = 4, size = 0.5) +
  #theme_bw()

6. Submit a link to this document in R Pubs to your Moodle. This assignment is worth 25 points.