Estimate the mean (expected value = \(E(Y_{ij})\)) of pefr based on the Wright peak flow meter rate with 95% CI.

Table 1 below shows the expected grand mean of PEFR for all groups.

cap <- "PEFR grand mean estimate along with 95% CI"
emmeans(m, ~ 1) %>% 
      tidy %>% 
      select(estimate, conf.low, conf.high) %>% 
      set_names(c("Estimate", "Upper", "Lower")) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 2))) %>% 
      kable_styling() %>% 
      add_header_above(c("", "95% Confidence Interval" = 2))

Estimate the variance \(var(Y_{ij})\) of pefr based on the Wright peak flow meter rate.

Table 2 shows the expected overall variance in measurements of PEFR

cap <- "Variance of PEFR"
emmeans(m, ~ 1) %>% 
      tidy %>% 
      select(std.error) %>% 
      set_names(c("Variance")) %>% 
      mutate(Variance = Variance^2) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = "c") %>% 
      kable_styling()

Estimate the covariance \(cov(Y_{ij},Y_{i'j})\) of pefr based on the Wright peak flow meter rate between two occasions.

Table 3 shows the expected covariance in measurements of PEFR on the same individual.

cap <- "Variance and covariance components of model"
tidy(m)[3:4, c(2, 4)] %>% 
      spread(group, estimate) %>% 
      select(Residual) %>% 
      mutate(Residual = Residual^2) %>% 
      set_names(c("$\\hat\\sigma^2_{id}$")) %>% 
      kable(cap = table_cap(cap),
            digits = 3,
            escape = F,
            align = "c") %>%
      kable_styling()

Estimate the correlation \(cor(Y_{ij},Y_{i'j})\) of pefr based on the Wright peak flow meter rate between two occasions. This correlation is called the intraclass correlation.

Table 4 shows the calculated intraclass correlation (ICC) between measurements of PEFR on the same individual.

cap <- "Intraclass correlation"
tidy(m)[3:4, c(2, 4)] %>% 
      spread(group, estimate) %>% 
      mutate(id = id^2,
             Residual = Residual^2,
             covar = id / (id + Residual)) %>% 
      set_names(c("$\\hat\\sigma^2_{id}$", "$\\hat\\sigma^2$", "ICC")) %>% 
      kable(cap = table_cap(cap),
            digits = 3,
            escape = F,
            align = c("l", rep("c", 2))) %>%
      kable_styling()

Create a “spaghetti” plot of the data.

d %>% 
      ggplot(aes(lrt, gcse, group = school)) +
      geom_line()

Figure 1: Observed trajectories for all schools

Does lrt predict the outcome of the gcse and does it depend on school type?

Table 5 below shows an F test for the interaction between LRT and School type which turned out not to be significant.

cap <- " "
m <- lmer(gcse ~ lrt + schgend + lrt : schgend + (lrt | school), data = d)
drop1(m) %>% 
      tidy %>% 
      select(term, statistic, p.value) %>% 
      set_names(c("Term", "F", "p")) %>% 
      kable(digits = 3,
            caption = table_cap(cap),
            align = c("l", rep("c", 2))) %>% 
      kable_styling()

The interaction effect was removed from the model and another round of tests was performed with only the marginal effects. Table 6 shows that both main effects should be in the model.

cap <- " "
m <- lmer(gcse ~ lrt + schgend + (lrt | school), data = d)
drop1(m) %>% 
      tidy %>% 
      select(term, statistic, p.value) %>% 
      mutate(term = fct_recode(term,
                               "LRT" = "lrt",
                               "School type" = "schgend")) %>% 
      set_names(c("Term", "F", "p")) %>% 
      kable(digits = 3,
            caption = table_cap(cap),
            align = c("l", rep("c", 2))) %>% 
      kable_styling()

Table 7 below shows fitted parameters for the mixed linear model.LRT and School type were fitted as fixed effects, and SchoolID was fitted as a random effect. A higher LRT score seems to indicate a higher value of GCSE. It also seems that all-girls schools have the highest GCSE after correcting for LRT

tidy(m) %>% 
      mutate(effect = fct_recode(effect,
                                 "Fixed" = "fixed",
                                 "Random" = "ran_pars"),
             group = fct_recode(group,
                                "School" = "school",
                                "Residual" = "Residual"),
             term = fct_recode(term, 
                               "Intercept" = "(Intercept)",
                               "LRT$^1$" = "lrt",
                               "School type: (Boys - Mixed)" = "schgend2",
                               "School type: (Girls - Mixed)" = "schgend3",
                               "LRT$^1$ $\\times$ (Boys - Mixed)" = "lrt:schgend2",
                               "LRT$^1$ $\\times$ (Girls - Mixed)" = "lrt:schgend3",
                               "$\\hat\\sigma:$ Intercept" = "sd__(Intercept)",
                               "$\\hat\\sigma:$ LRT$^1$" = "sd__lrt",
                               "$\\hat\\rho:$ Intercept $\\times$ LRT$^1$" = "cor__(Intercept).lrt",
                               "$\\hat\\sigma:$ Residual" = "sd__Observation")) %>%
      set_names(c("Effect", "Group", "Term", "Estimate", "SE", "t", "DF", "p")) %>%
      arrange(Effect, Group, Term) %>% 
      select(-Effect, -Group) %>% 
      kable(caption = table_cap("Summary table for model parameters"),
            digits = 3,
            escape = FALSE,
            align = c("l", rep("c", 5))) %>% 
      kable_styling(bootstrap_options = c("striped", "hover")) %>% 
      group_rows(index = c("", "London Reading Test" = 1, "School Type" = 2, "Random Effects" = 4)) %>% 
      footnote(number = c("Normalized by subtracting the mean and dividing by the standard deviation"))

Is there a signfificant variability on how lrt predict gcse between schools?

Table 8 shows the outcome of comparing models with and without a random effect for the variability in lrt between schools. As can be seen the test score is high and \(p < 0.0005\) indicating that there is variability between schools.

cap <- " "
ranova(m) %>% 
      tidy %>% 
      select(-npar, -logLik, -df) %>% 
      set_names(c("Term", "AIC", "LRT", "p")) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 2))) %>% 
      kable_styling()

Write out your final model using greek letters for parameters associated with fixed effects and small cap roman letters for random coefficients.

\[ \begin{aligned} Y_{GCSE} = -&0.959 \\ +&5.581 \cdot LRT \\ +&0.934 \cdot \mathbb I_{\text{All-Boys}} \\ +&2.598 \cdot \mathbb I_{\text{All-Girls}} \\ +&s_{\text{Intercept}} \\ +&s_{\text{LRT}} \\ +&s_{\text{Residual}} \\s_{\text{Intercept}} = \quad &\mathcal N(0, 2.816^2) \\s_{\text{LRT}} = \quad &\mathcal N(0, a^2) \\a = \quad &1.190 \cdot LRT \\s_{\text{Residual}} = \quad &\mathcal N(0, 7.443^2) \end{aligned} \]

Based on your model: Find the 5 schools with highest and lowest gcse score adjusted to lrt=0 for each shcool type.

Table 9 below shows the bottom five schools for each type and table 10 below that shows the top five.

cap <- "Bottom 5 schools for each type"
d %>% 
      distinct(school, schgend) %>% 
      mutate(lrt = 0) %>% 
      add_predictions(m) %>% 
      group_by(schgend) %>% 
      top_n(-5, pred) %>% 
      select(schgend, school, pred) %>% 
      arrange(schgend, pred) %>% 
      set_names(c("School Type", "School", "Prediction")) %>% 
      kable(cap = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 2))) %>% 
      kable_styling(bootstrap_options = c("striped", "hover")) %>% 
      collapse_rows(1)

cap <- "Top 5 schools for each type"
d %>% 
      distinct(school, schgend) %>% 
      mutate(lrt = 0) %>% 
      add_predictions(m) %>% 
      group_by(schgend) %>% 
      top_n(5, pred) %>% 
      select(schgend, school, pred) %>% 
      arrange(schgend, desc(pred)) %>% 
      set_names(c("School Type", "School", "Prediction")) %>% 
      kable(cap = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 2))) %>% 
      kable_styling(bootstrap_options = c("striped", "hover")) %>% 
      collapse_rows(1)

Describe the structure of the data.

The data are composed of 6 repeated measurements taken from each of 14 different loblolly pine trees at 5-year increments.

cap <- "View of data structure."
d %>% 
      select(seed, age, height) %>% 
      spread(age, height) %>% 
      rename("Tree Seed" = seed) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 6))) %>% 
      kable_styling(bootstrap_options = c("hover", "striped")) %>% 
      add_header_above(c("", "Age (in decades)" = 6))

Provide mean height by time and variance by time. Comment on how the mean and variance change with time.

Figure 2 below shows the mean and variance of tree heights over time. The mean height seems to be mostly linear but with a small negative second-order term. Seeing such a slow decrease in growth could stem from trees continuing to grow for many years and the sample is recorded for only 25 years. The variance seems to grow linearly with age with a discontinuity around 20 years. The reduction in variance might stem from a decreased growth-rate. As the trees grow more slowly the difference between them should also grow more slowly.

d %>% 
      group_by(age) %>% 
      summarize(Mean = mean(height),
                Variance = var(height)) %>% 
      gather(variable, value, -age) %>% 
      ggplot(aes(age, value)) +
      geom_line() +
      facet_wrap("variable", scales = "free") +
      labs(x = "Age",
           y = "Value")

Figure 2: Mean and variance of height over time.

Provide a “spaghetti plot” of the data.

Figure fig_num + 1 below shows a spaghetti plot of the observed growth trajectories.

d %>% 
      ggplot(aes(age, height, group = seed)) +
      geom_line()

Fit the appropriate model using polynomial regression (no higher than third degree).

We start out with a complicated model containing fixed effects up to the third degree and second degree random effects.

m <- lmer(height ~ age_cent + I(age_cent^2) + I(age_cent^3) + (age_cent + I(age_cent^2) | seed), data = d)

Simplifying variance structure

We first simplify the variance structure as needed. Table table_num + 1 below shows that no random effects should be removed.

cap <- "Likelihood ratio tests for random effects"
ranova(m) %>% 
      tidy %>% 
      select(-npar, -logLik, -df) %>% 
      set_names(c("Term", "AIC", "LRT", "p")) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 3))) %>% 
      kable_styling()

cap <- "Significance tests of fixed effects"
Anova(m, type = "III", test = "F") %>% 
      tidy %>% 
      select(-df, -Df.res) %>% 
      set_names(c("Term", "F Statistic", "p")) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 2))) %>% 
      kable_styling()

m <- lmer(height ~ age_cent + I(age_cent^2) + (age_cent + I(age_cent^2) | seed), data = d)

cap <- "Significance tests for fixed effects"
Anova(m, type = "III", test = "F") %>% 
      tidy %>% 
      select(-df, -Df.res) %>% 
      set_names(c("Term", "F Statistic", "p")) %>% 
      kable(caption = table_cap(cap),
            digits = 3,
            align = c("l", rep("c", 2))) %>% 
      kable_styling()

tidy(m) %>% 
      mutate(effect = fct_recode(effect,
                                 "Fixed" = "fixed",
                                 "Random" = "ran_pars"),
             group = fct_recode(group,
                                "School" = "school",
                                "Residual" = "Residual"),
             term = fct_recode(term, 
                               "Intercept$^a$" = "(Intercept)",
                               "Age (In decades)$^b$" = "age_cent",
                               "Age$^2$" = "I(age_cent^2)",
                               "$\\hat\\sigma:$ Intercept" = "sd__(Intercept)",
                               "$\\hat\\sigma:$ Age" = "sd__age_cent",
                               "$\\hat\\sigma:$ Age$^2$" = "sd__I(age_cent^2)",
                               "$\\hat\\rho:$ Intercept $\\times$ Age" = "cor__(Intercept).age_cent",
                               "$\\hat\\rho:$ Intercept $\\times$ Age$^2$" = "cor__(Intercept).I(age_cent^2)",
                               "$\\hat\\rho:$ Age $\\times$ Age$^2$" = "cor__age_cent.I(age_cent^2)",
                               "$\\hat\\sigma:$ Residual" = "sd__Observation")) %>%
      set_names(c("Effect", "Group", "Term", "Estimate", "SE", "t", "DF", "p")) %>%
      select(-Effect, -Group) %>% 
      kable(caption = table_cap("Summary table for model parameters"),
            digits = 3,
            escape = FALSE,
            align = c("l", rep("c", 5))) %>% 
      kable_styling(bootstrap_options = c("striped", "hover")) %>% 
      group_rows(index = c("", "Age" = 2, "Random Effects" = 7)) %>% 
      footnote(alphabet = c("Height at age = 13 years", 
                            "Centered by subtracting its mean (13 years) before modelling."))

Comment on the shape of the growth curve of the final model and derive the formula for growth.

The estimated growth curve has a positive linear shape with a negative second-order term, implying that the growth slows down over time.

\[ \begin{aligned} x &= Age - 1.3 \\ Y_{\text{Height}} &= 35.438 + x \cdot \beta_\text{Growth}\\ \beta_\text{Growth} &= 26.633 - 4.985 \cdot x + s^2_{\text{Growth}} \\ s^2_{\text{Growth}} &= \mathcal N(0, a^2) \\ a^2 &= \sigma^2_{\text{Seed}} + \sigma^2_{\text{Residual}} \\ & \hspace{3.4em} + x (x \cdot \sigma^2_{x} + 2\sigma_{\text{Seed}\times x}) \\ &\hspace{3.4em}+ x^2 (x^2 \cdot \sigma^2_{x^2} + 2\sigma_{\text{Seed}\times x^2}) \\ &\hspace{3.4em} + 2x^3 \cdot \sigma_{x\times x^2} \\ &= 3.07 + 0.350 \\ & \hspace{3.4em} + x (x \cdot 0.925+ 2 \cdot 1.51) \\ &\hspace{3.4em}+ x^2 (x^2 \cdot 0.382 - 2 \cdot 0.949) \\ &\hspace{3.4em} - x^3 \cdot 2 \cdot 0.342 \end{aligned} \]

Comment on the variance of the growth over time. In other words, what does the model you choose imply about the variance over time. In addition you could compare that to the sample estimates above.

Figure 4 shows the observed variance by age and the fitted variance estimated via the model’s covariance matrix.

var_func <- function(age) 3.07 + 0.35 + age * (age * 0.975 + 2 * 1.51) + age^2 * (age^2 * 0.382 - 2 * 0.949) - age^3 * 2 * 0.342
d %>% 
      group_by(age) %>% 
      summarize(Observed = var(height)) %>% 
      mutate(Fitted = var_func(age - 1.3)) %>% 
      gather(Type, value, -age) %>% 
      ggplot(aes(age, value, col = Type)) +
      geom_line() + 
      scale_color_jama(guide = guide_legend(title.position = "top", title.hjust = 0.5)) +
      labs(x = "Age",
           y = "Variance") +
      theme(legend.position = c(0.1, 0.8))

Figure 4: Observed and fitted variance by age

Provide the fitted mean growth vs. age and plot it with 95% confidence bands.

Figure 5 below shows the estimated mean growth along with 95% confidence bands.

emmeans(m,  ~ age_cent, cov.reduce = F) %>% 
      tidy(level = 0.95) %>% 
      ggplot(aes(age_cent + 1.3, estimate,
                 ymin = conf.low, 
                 ymax = conf.high)) + 
      geom_line() +
      geom_ribbon(alpha = 0.2) +
      labs(x = "Age",
           y = "Height")

Figure 5: Estimated mean growth and 95% confidence bands

Estimate the average growth the first decade, based on the final model you choose.

Table 16 shows the estimated growth in the first decade along with its 95% confidence interval.

emmeans(m, ~ age_cent, at = list(age = -0.3)) %>% 
      tidy %>% 
      select(estimate, std.error, conf.low, conf.high) %>% 
      set_names(c("Estimate", "SE", "Lower", "Upper")) %>% 
      kable(caption = table_cap("Estimated growth in the first decade"),
            digits = 3,
            align = c("l", "c", "c", "c")) %>% 
      kable_styling() %>% 
      add_header_above(c("", "", "95% Confidence Interval" = 2))

Look at residuals

Figure fig_num + 1 below shows the distribution of model residuals. They seem to be approximately normally distributed with mean zero. Figure fig_num + 2 shows a pairplot for estimated random effects. We can see the effects are approximately randomly distributed.

augment(m) %>% 
      ggplot(aes(.resid)) +
      geom_histogram(aes(y = ..density..)) +
      geom_density(fill = "grey50", alpha = 0.5) +
      labs(x = "Residual", y = "Density")

Figure 6: Distribution of residuals

ranef(m)$seed %>% 
      as_tibble %>% 
      ggpairs()

Figure 7: Pairplot of random effects

Select 12 trees randomly, and plot observed and fitted growth at the individual level using a facet grid.

Figure 8 shows observed and fitted growth for twelve randomly selected trees.

d %>% 
      filter(seed %in% sample(unique(seed), 12)) %>% 
      add_predictions(m) %>% 
      ggplot(aes(age, height, group = seed)) +
      geom_line(aes(age, pred)) +
      geom_point() +
      facet_wrap("seed") +
      labs(x = "Age", y = "Height")

Figure 8: Fitted growth curves along with observed trajectory values for twelve randomly sampled seeds.

The change in TOTALFLOW in each group with 95% CI and if it was statistically significant.

Table 17 shows the estimated contrast along with significance tests and 95% confidence intervals for the change in TOTALFLOW for each group.

emmeans(m, pairwise ~ visitf | af_statusf)$contrasts %>% 
      tidy %>% 
      unite(Contrast, level1, level2) %>% 
      mutate(Contrast = "Post - Pre",
             estimate = -estimate,
             upper = round(estimate + qt(0.975, df) * std.error, 3),
             lower = round(estimate - qt(0.975, df) * std.error, 3)) %>% 
      select(Contrast, af_statusf, estimate, lower, upper, p.value) %>% 
      set_names(c("Contrast", "Group", "Estimate", "Lower", "Upper", "p")) %>% 
      kable(caption = table_cap("Estimated pre-post change for each group"),
            digits = 3,
            align = c("l", rep("c", 5))) %>% 
      kable_styling() %>% 
      add_header_above(c("", "", "", "95% Confidence Interval" = 2, ""))

The estimated difference in the change between groups with 95%, and if the change in the SR group was significantly different from the AF group. Use emmeans.

Table 18 shows the estimated difference between groups in TOTALFLOW change along with significance test and 95% confidence interval.

tidy(m, effects = "fixed", conf.int = TRUE) %>% 
      filter(str_detect(term, ":")) %>% 
      select(term, estimate, conf.low, conf.high, p.value) %>% 
      mutate(term = "AF - SR") %>% 
      set_names(c("Contrast", "Estimate", "Lower", "Upper", "p")) %>% 
      kable(caption = table_cap("Estimated difference in pre-post change between groups"),
            digits = 3,
            align = c("l", rep("c", 4))) %>% 
      kable_styling() %>% 
      add_header_above(c("", "", "95% Confidence Interval" = 2, ""))

A short description of the methodology for his paper.

In this clinical trial participant levels of brain flow were measured before and after a treatment. Participants were grouped by whether the treatment was deemed successful (SR group) or not (AF group). A linear mixed model (LMM) was used to estimate the treatment effect as well as any difference in effect between the two groups. The LMM was fitted according to the formula

\[ Y_{ij} = \beta_0 + \beta_1 X_i + \beta_2 t_{ij} + \beta_3 t_{ij} X_i + \beta_4S_i + \beta_5A_i + \epsilon_{ij}, \]

where \(Y_{ij}\) is the j’th measurement for subject i, \(X_i\) is the treatment group of subject i, \(t_ij\) indicates whether \(Y_{ij}\) is a pre- or post measurement for subject i, and \(S_i\) and \(A_i\) are the sex and age of subject i. The estimated difference in pre-post change between groups is then modeled by the parameter \(\beta_3\). Age and sex were included in the model to adjust for any potantial confounding.

The spaghetti plot

Figure 9 shows the observed trajectories for participants grouped by whether their treatment was deemed successful or unsuccessful.

d %>% 
      select(id, totalflow, visitf, af_statusf) %>% 
      mutate(visitf = fct_recode(visitf,
                                 "Pre" = "1",
                                 "Post" = "2")) %>% 
      ggplot(aes(visitf, totalflow, group = id)) +
      geom_line() +
      facet_grid("af_statusf")

Figure 9: Observed trajectories for participants

A look at total residuals.

Figure fig_num + 1 below shows the distribution of model residuals. They seem to be approximately normally distributed with mean zero. Figure fig_num + 2 shows the distribution of estimated random intercepts. We can see the intercepts are approximately randomly distributed.

augment(m) %>% 
      ggplot(aes(.resid)) +
      geom_histogram(aes(y = ..density..)) +
      geom_density(fill = "grey50", alpha = 0.5) +
      labs(x = "Residual", y = "Density")

Figure 10: Distribution of residuals

ranef(m)[[1]] %>% 
      as_tibble %>% 
      set_names(c("intercept")) %>% 
      ggplot(aes(intercept)) +
      geom_histogram(aes(y = ..density..)) +
      geom_density(fill = "grey50", alpha = 0.5) +
      labs(x = "Intercept", y = "Density")

Figure 11: Distribution of random intercepts

Compute the contrast in part 2) using glht of multcomp and comment on the difference.

The main difference between using emmeans and glht is that glht requires the use of matrices to calculate your linear hypothesis.

glht(m, matrix(c(0, 0, 0, 0, 0, 1), ncol = 6, byrow = T)) %>% 
      summary() %>% 
      tidy %>% 
      mutate(contrast = "AF - SF") %>% 
      select(contrast, estimate, std.error, p.value) %>%
      set_names(c("Contrast", "Estimate", "SE", "p")) %>% 
      kable(digits = 3,
            caption = table_cap("Estimated difference in pre-post change between groups using glht")) %>% 
      kable_styling()