Props_mixed_model.Rmd

---
title: "Mixed effect model analyses"
output: html_notebook
---

```{r setup_Feates_mixedmod, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Mixed model analysis of electrophysiological features


Core analyses used for investigation of inter-animal variability in intrinsic properties of layer 2 stellate cells and their dorsoventral organisation. This code uses functions in "Functions.Rmd".

Organisation:

Fit mixed models
Compare fits of mixed models to one another and to population level linear model
Extract other summary data from the models
Generate summary table
Plot fits of mixed models
Generate plots to evaluate assumptions of mixed models



## Fit mixed models

Fit mixed models to all measured properties using lmer. Models are fit with a random intercept and slope (_vsris), a random intercept only (_vsri) and with a correlated random intercept and slope (_vscris). Each model is described in its corresponding function, e.g. model_vsris, and has a corresponding null model, e.g. model_vsris_null. Fits with _lT use lmerTest to obtain alternative significance estimates.
```{r}
data.sc_r <- data.sc_r %>%
  mutate(mm_vsris = map(data, model_vsris))%>%
  mutate(mm_vsris_null = map(data, model_vsris_null))%>%
  mutate(mm_vsri = map(data, model_vsri))%>%
  mutate(mm_vsri_null = map(data, model_vsri_null))%>%
  mutate(mm_vscris = map(data, model_vscris))%>%
  mutate(mm_vscris_null = map(data, model_vscris)) %>%
  mutate(mm_vsris_lT = map(data, model_vsris_lT))
```
Similar to above, but additional models using the transformed properties, or the transformed properties and transformed dvlocmm (vsris_1). Transformations were carried out in LoadData.Rmd. 
```{r}
data.sc_TP_r <- data.sc_TP_r %>%
  mutate(mm_vsris = map(data, model_vsris))%>%
  mutate(mm_vsris_null = map(data, model_vsris_null))%>%
  mutate(mm_vsris_1 = map(data, model_vsris_1))%>%
  mutate(mm_vsris_1_null = map(data, model_vsris_1_null))
```




## Compare fits of mixed models to one another and to population level linear model

Extract AIC for all models.
```{r Extract AIC}
data.sc_r <- data.sc_r %>% 
  mutate(vsris_glance = map(mm_vsris, broom::glance)) %>%
  mutate(vsris_null_glance = map(mm_vsris_null, broom::glance)) %>%
  mutate(vsri_glance = map(mm_vsri, broom::glance)) %>%
  mutate(vscris_glance = map(mm_vscris, broom::glance)) %>%
  mutate(vsris_AIC = map_dbl(vsris_glance, ~.$AIC)) %>%
  mutate(vsris_null_AIC = map_dbl(vsris_null_glance, ~.$AIC)) %>%
  mutate(vsri_AIC = map_dbl(vsri_glance, ~.$AIC)) %>%
  mutate(vcsris_AIC = map_dbl(vscris_glance, ~.$AIC))

data.sc_TP_r <- data.sc_TP_r %>%
  mutate(vsris_glance = map(mm_vsris, broom::glance)) %>%
  mutate(vsris_null_glance = map(mm_vsris_null, broom::glance)) %>%
  mutate(vsris_1_glance = map(mm_vsris_1, broom::glance)) %>%
  mutate(vsris_1_null_glance = map(mm_vsris_1_null, broom::glance)) %>%
  mutate(vsris_AIC = map_dbl(vsris_glance, ~.$AIC)) %>%
  mutate(vsris_null_AIC = map_dbl(vsris_null_glance, ~.$AIC)) %>%
  mutate(vsris_1_AIC = map_dbl(vsris_1_glance, ~.$AIC)) %>%
  mutate(vsris_1_null_AIC = map_dbl(vsris_1_null_glance, ~.$AIC))
```

In general vsris and vscris have similar AIC (compare columns in data.sc_r). We focus the subsequent analysis on vsris. For discussion of maximal models versus random intercept only models see Barr et al. Journal of Memory and Language, 2013.


## Example data suggests substantial inter-animal variability

Plot example animals
```{r}
(rheo_example <- ggplot(filter(data.sc, id == "mouse_20131106" | id == "mouse_20140113"), aes(dvloc, rheo)) +
  geom_point(aes(colour = id), shape = 3) +
  labs(x = "Location (µm)", y = "Rheobase (pA)", colour = "Mouse")) +
  scale_colour_manual(labels = c("20131106", "20140113"), values = c("red", "blue"))
(resf_example <- ggplot(filter(data.sc, id == "mouse_20131106" | id == "mouse_20140113"), aes(dvloc, resf)) +
  geom_point(aes(colour = id)) +
  labs(x = "Location (µm)", y = "Resonance frequency (Hz)") +
  theme(legend.position = ""))
```




Make plots for fits using linear models for each id or mixed models. First look at input resistance.
```{r}
ir_predictplot <- ggplot(filter(data.sc_r, property == "ir")$data[[1]], aes(x = dvlocmm, y = value, group = id)) +
    stat_smooth(geom = "line", method = lm, se = FALSE)
ir_predictplot <- gg_ir_format(ir_predictplot, 0, 60)

mm_ir_predictplot <- predict_plot(filter(data.sc_r, property == "ir")$data[[1]],
                                  filter(data.sc_r, property == "ir")$mm_vsris[[1]])
mm_ir_predictplot <- gg_ir_format(mm_ir_predictplot, 0, 60)

ir_predictplot
mm_ir_predictplot
```

Next look at rheobase.
```{r}
rheo_predictplot <-  ggplot(filter(data.sc_r, property == "rheo")$data[[1]], aes(x = dvlocmm, y = value, group = id)) +
    stat_smooth(geom = "line", method = lm, se = FALSE)
rheo_predictplot <- gg_rheo_format(rheo_predictplot, 0, 600)

mm_rheo_predictplot <- predict_plot(filter(data.sc_r, property == "rheo")$data[[1]],
                               filter(data.sc_r, property == "rheo")$mm_vsris[[1]])
mm_rheo_predictplot  <- gg_rheo_format(mm_rheo_predictplot, 0, 600)

rheo_predictplot
mm_rheo_predictplot

```


And then look at resonant frequency.
```{r}
resf_predictplot <-  ggplot(filter(data.sc_r, property == "resf")$data[[1]], aes(x = dvlocmm, y = value, group = id)) +
    stat_smooth(geom = "line", method = lm, se = FALSE)
resf_predictplot <- gg_resf_format(resf_predictplot, 0, 12)


mm_resf_predictplot <- predict_plot(filter(data.sc_r, property == "resf")$data[[1]],
                               filter(data.sc_r, property == "resf")$mm_vsris[[1]])

mm_resf_predictplot  <- gg_resf_format(mm_resf_predictplot, 0, 12)

resf_predictplot
mm_resf_predictplot
```





## Extract summary statistics from models.

The function mixedmod_extract is used to return statistics for mm_vsris. Focus on the model with random intercept and slope (mm_vsris). Store model gradient (extracted with summary / glance), marginal and conditional R2 (extracted with r.squaredGLMM) and p-value vs null model (calculated with ANOVA vs null model). Also extract model slopes. Creation of data.sc_r_1 is for separate handling of results with transformed dvlocmm.
```{r Extract Summary statistics, warning=FALSE}
data.sc_r <- mixedmod_extract(data.sc_r, mm_vsris)
data.sc_TP_r <- mixedmod_extract(data.sc_TP_r, mm_vsris)
data.sc_TP_r <- mixedmod_extract(data.sc_TP_r, mm_vsris_1)
```


## Compare fits of mixed models to population level linear model

To test whether effects of animal id are significant compare mixed model fits with linear model fits. Modified from: https://web.stanford.edu/class/psych252/section/Mixed_models_tutorial.html.
```{r Compare mixed with linear model using chisq}
## linearmodel_to_fit fits: lm(value ~ dvlocmm, data = df, na.action = na.exclude)
data.sc_r <- data.sc_r %>%
  mutate(linearmodel = map(data, linearmodel_to_fit))

data.sc_r <- mixed_vs_linear_pchisqu(data.sc_r, mm_vsris, linearmodel) 
data.sc_r$mm_vsris_vslinear_pdiff_adj <- p.adjust(data.sc_r$mm_vsris_vslinear_pdiff, method = "BH")

# Look at linear model output
select(data.sc_r, property, linearmodel) %>%
  mutate(glance = map(linearmodel, broom::glance),
         adj.r2 = map(glance, ~.$adj.r.squared)) %>%
  select(property, adj.r2) %>%
  kableExtra::kable()

table_mixedvslinear(data.sc_r, "mm_vsris_vslinear", "property")
```

And, do the same thing for the transformed data:
```{r}
data.sc_TP_r <- data.sc_TP_r %>%
  mutate(linearmodel = map(data, linearmodel_to_fit)) %>%
  mutate(linearmodel_1 = map(data, linearmodel_to_fit_1))

data.sc_TP_r <- mixed_vs_linear_pchisqu(data.sc_TP_r, mm_vsris, linearmodel)
data.sc_TP_r$mm_vsris_vslinear_pdiff_adj <- p.adjust(data.sc_TP_r$mm_vsris_vslinear_pdiff, method = "BH")

table_mixedvslinear(data.sc_TP_r, "mm_vsris_vslinear", "property")

data.sc_TP_r <- mixed_vs_linear_pchisqu(data.sc_TP_r, mm_vsris_1, linearmodel_1)
data.sc_TP_r$mm_vsris_1_vslinear_pdiff_adj <- p.adjust(data.sc_TP_r$mm_vsris_1_vslinear_pdiff, method = "BH")

table_mixedvslinear(data.sc_TP_r, "mm_vsris_1_vslinear", "property")
```

Alternative approaches to compare the mixed model with the reduced model. See: http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-significance-of-random-effects. Example below uses ANOVA directly. Note, this doesn't work when calling the model stored within data.sc_r (commented out).
```{r, warning=FALSE}
#a <- data.sc_r$linearmodel[[1]]
#b <- data.sc_r$mm_vsris[[1]]
a <- lmer(vm ~ dvlocmm + (dvlocmm||id), data = data.sc)
b <- lm(vm ~ dvlocmm, data = data.sc)
anova(a, b)
```



## Compare the model with a null model containing only id as a random effect
Use ANOVA to compare the model with the null model.
To do: compare with car::Anova.
```{r Compare to null model, warning=FALSE}
data.sc_r <- data.sc_r %>%
  mutate(anova = map2(mm_vsris, mm_vsris_null, ~anova(.x,.y))) %>%
  mutate(tidy_anova = map(anova, broom::tidy)) %>%
  mutate(anova_p_val = map_dbl(tidy_anova, ~.$p.value[2]))

data.sc_r <- data.sc_r %>%
  mutate(anova_p_val_adj = p.adjust(data.sc_r$anova_p_val, method = "BH"))
```
And, again for the transformed data:
```{r, warning=FALSE}

data.sc_TP_r <- data.sc_TP_r %>%
  mutate(anova = map2(mm_vsris, mm_vsris_null, ~anova(.x,.y))) %>%
  mutate(tidy_anova = map(anova, broom::tidy)) %>% 
  mutate(anova_p_val = map_dbl(tidy_anova, ~.$p.value[2]))
data.sc_TP_r <- data.sc_TP_r %>%
  mutate(anova_p_val_adj = p.adjust(data.sc_r$anova_p_val, method = "BH"))

data.sc_TP_r <- data.sc_TP_r %>%
  mutate(anova_1 = map2(mm_vsris_1, mm_vsris_1_null, ~anova(.x,.y))) %>%
  mutate(tidy_anova_1 = map(anova_1, broom::tidy)) %>%
  mutate(anova_p_val_1 = map_dbl(tidy_anova_1, ~.$p.value[2]))
data.sc_TP_r <- data.sc_TP_r %>%
  mutate(anova_p_val_adj_1 = p.adjust(data.sc_TP_r$anova_p_val_1, method = "BH"))
```


## Generate summary tables

props_table_unnest

Parameters from fitting models directly to the experimentally obtained properties.
```{r Table_mixed_model_properties}
props_for_table <- c("property", "anova_p_val", "anova_p_val_adj", "mm_vsris_marginal.r2", "mm_vsris_conditional.r2", "mm_vsris_vslinear_pdiff", "mm_vsris_vslinear_pdiff_adj", "mm_vsris_gradient_slopes", "mm_vsris_slope_min", "mm_vsris_slope_max")
props_table_unnest <- unnest(data.sc_r[props_for_table])

props_table_unnest %>%
  knitr::kable(
    digits = 5,
    caption = "Fit of measured membrane properties as a function of location",
    col.names = c("property", "p (ANOVA)", "p_adj (ANOVA)", "marginal R2", "conditional R2", "p (vs linear)", "p_adj (vs linear)", "slope", "slope (min)", "slope (max)")
) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))

write_csv(props_table_unnest, "results_model_table.csv")

# reduced version for the manuscript
props_man_table <- c("property", "mm_vsris_gradient_slopes", "anova_p_val_adj", "mm_vsris_marginal.r2", "mm_vsris_conditional.r2", "mm_vsris_slope_min", "mm_vsris_slope_max", "mm_vsris_vslinear_pdiff_adj")
man_table <- unnest(data.sc_r[props_man_table])

man_table$property <- c("Vm (mV)", "IR (MΩ)", "Sag", "Tm (ms)", "Res. frequency (Hz)", "Res. magnitude", "Spike thresold (mV)", "Spike maximum (mV)", "Spike width (ms)", "Rheobase (pA)", "Spike AHP (mV)", "I-F slope (Hz/pA)")

man_table <- man_table %>%
  rownames_to_column() %>%
  arrange(desc(mm_vsris_marginal.r2)) %>%
  column_to_rownames() %>%
  as_tibble()

man_table$anova_p_val_adj <- format(man_table$anova_p_val_adj, digits = 3)
man_table$mm_vsris_vslinear_pdiff_adj <- format(man_table$mm_vsris_vslinear_pdiff_adj, digits = 3)

(man_table <- man_table %>%
  knitr::kable(
    digits = 3,
    col.names = c("Feature", "Slope", "p (slope)", "Marginal R2", "Conditional R2", "Slope (min)", "Slope (max)", "p (vs linear)")
) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover")))


if (stajpeg == 1) {
  table_path <- paste0(getwd(), "/Tables/mm_prop_vs_dv_summary.jpg")
   man_table %>% kableExtra::save_kable(file = table_path, self_contained = T)
}

table_path <- paste0(getwd(), "/Tables/mm_prop_vs_dv_summary.html")
man_table %>% kableExtra::save_kable(file = table_path, self_contained = T)
```

Parameters from fitting models with properties transformed.
```{r Table_mixed_model_properties_t}
props_table_t <- as.tibble(data.sc_TP_r[props_for_table])
props_table_t_unnest <- unnest(props_table_t)

props_table_t_unnest %>%
  knitr::kable(
  digits = 5,
  caption = "Fit of measured membrane properties as a function of location",
    col.names = c("property", "p (ANOVA)", "p_adj (ANOVA)", "marginal R2", "conditional R2", "p (vs linear)", "p_adj (vs linear)", "slope", "slope (min)", "slope (max)")
) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
```

Parameters from fitting models with properties and dvlocmm transformed.
```{r Table_mixed_model_properties_1}
props_for_table <- c("property", "anova_p_val_1", "anova_p_val_adj_1", "mm_vsris_1_marginal.r2", "mm_vsris_1_conditional.r2", "mm_vsris_1_vslinear_pdiff", "mm_vsris_1_vslinear_pdiff_adj", "mm_vsris_1_gradient_slopes", "mm_vsris_1_slope_min", "mm_vsris_1_slope_max")

props_table_1 <- as.tibble(data.sc_TP_r[props_for_table])
props_table_1_unnest <- unnest(props_table_1)

props_table_1_unnest %>%
  knitr::kable(
  digits = 5,
  caption = "Fit of measured membrane properties as a function of location",
    col.names = c("property", "p (ANOVA)", "p_adj (ANOVA)", "marginal R2", "conditional R2", "p (vs linear)", "p_adj (vs linear)", "slope", "slope (min)", "slope (max)")
) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
```



## Plot fits of mixed models

We want to plot for each model the prediction at location = 0 for each animal (I), the model prediction for location = 1 mm (I + S) and a line indicating the slope with start centred at the value of the population level model at location = 0.

Reformat the data to generate plots.
Call to prep_int_slopes extracts model predictions ready for plotting.
```{r Format mixed model outputs ready for plotting, warning=FALSE}
combined_intercepts_slopes <- prep_int_slopes(data.sc_r, "property", "mm_vsris")

id_housing <-  distinct(data.sc, id, housing, age, sex)
combined_intercepts_slopes <- left_join(combined_intercepts_slopes, id_housing, by = "id")

combined_intercepts_slopes$property_factors <- as.factor(combined_intercepts_slopes$property)

combined_intercepts_slopes$property_factors = factor(combined_intercepts_slopes$property_factors, c("vm", "ir", "sag", "tau", "resf", "resmag", "rheo", "fi", "ahp", "spkmax", "spkthr", "spkhlf"))
```


Now generate the plot.
```{r Facetted_plot_of_model_fits}
labels_intercepts <- c(ahp = "AHP min. (mV)", fi = "F-I (Hz / pA)", ir = "IR (MΩ)", resf = "Res F (Hz)", resmag = "Res. mag.", rheo = "Rheobase (pA)", sag = "Sag", spkhlf = "Spike h-w (ms)", spkmax = "Spike max. (mV)", spkthr = "Spike thres. (mV)", tau = "Tm (ms)", vm = "Vrest (mV)")

IS_figure <- ggplot(combined_intercepts_slopes, aes(x = measure, y = value_1, colour = housing)) +
  geom_line(aes(group = id)) +
  geom_jitter(aes(y = value_2), width = 0.2, alpha = 0.5) +
  #geom_boxplot(aes(y = value_2), alpha = 0.1) +
  stat_summary(aes(y = value_2), fun.y = "mean", fun.ymin = "mean", fun.ymax = "mean", size = 0.2, geom = "crossbar") +
  scale_x_discrete(limits = c("ind_intercept", "ind_intercept_slope", "global_intercept", "global_intercept_slope"), label = c("I", "I + S", "", "")) +
  facet_wrap(~property_factors, scales = "free",  labeller = labeller(property_factors = labels_intercepts)) +
  theme_classic() +
  hist_theme +
  theme(axis.line.x = element_blank(), axis.ticks.x = element_blank(), legend.position = "bottom")

IS_figure
```



Make combined figure.
```{r warning=FALSE}
ggdraw() +
  draw_plot(rheo_predictplot, x = 0, y = .5, width = .18, height = .4) +
  draw_plot(mm_rheo_predictplot, x = 0, y = 0, width = .18, height = .4) +
  draw_plot(resf_predictplot, x = 0.2, y = .5, width = .18, height = .4) +
  draw_plot(mm_resf_predictplot, x = 0.2, y = 0, width = .18, height = .4) +
  draw_plot(IS_figure, x = 0.4, y = 0, width = .6, height = 1) +
  draw_plot_label(label = c("A", "B", "C"), size = 15,
                  x = c(0, 0, 0.4), y = c(1, 0.5, 1)) +
  draw_plot_label(label = c("Independent fits", "Mixed model fits"), size = 12,
                  x = c(0.006, 0.006), y = c(0.99, 0.49))
  
  
```



Save the figure.
```{r Save plot of model fits}
ggsave("Figures/I_S_figure.png", width = 220, height = 120, units = "mm")
```



## Generate plots to evaluate assumptions of mixed models
To test for linearity, and to assess homoscedasticity, plot residuals (.residual) versus fitted values (.fitted) generated by broom::augment.

```{r}
resid_plot_data <- unnest(data.sc_r, mm_vsris_aug) %>%
  select(property, id, .resid, .fitted)

ggplot(resid_plot_data, aes(.fitted, .resid)) +
  geom_point() +
  facet_wrap(~property, scales = "free") +
  theme_classic()
```

Make the same plot, but for transformed data
```{r}
resid_plot_data_TP <- unnest(data.sc_TP_r, mm_vsris_aug) %>%
  select(property, id, .resid, .fitted)

ggplot(resid_plot_data_TP, aes(.fitted, .resid)) +
  geom_point() +
  facet_wrap(~property, scales = "free") +
  theme_classic()
```

And for transformed location and data.
```{r}
resid_plot_data_TP_1 <- unnest(data.sc_TP_r, mm_vsris_1_aug) %>%
  select(property, id, .resid, .fitted)

ggplot(resid_plot_data_TP_1, aes(.fitted, .resid)) +
  geom_point() +
  facet_wrap(~property, scales = "free") +
  theme_classic()
```



Evaluate normality of residuals.
```{r}
ggplot(resid_plot_data, aes(.resid)) +
  geom_histogram() +
  facet_wrap(~property, scales = "free") +
  theme_classic()
```

Make same plot but for transformed data.
```{r}
ggplot(resid_plot_data_TP, aes(.resid)) +
  geom_histogram() +
  facet_wrap(~property, scales = "free") +
  theme_classic()
```
And for transformed location and data.
```{r}
ggplot(resid_plot_data_TP_1, aes(.resid)) +
  geom_histogram() +
  facet_wrap(~property, scales = "free") +
  theme_classic()
```




Q-Q plots of residuals for raw data
```{r}
ggplot(resid_plot_data, aes(sample = .resid)) + stat_qq() + stat_qq_line() + facet_wrap(~property, scales = "free")
```

Make same plots but for transformed data
```{r}
ggplot(resid_plot_data_TP, aes(sample = .resid)) + stat_qq() + stat_qq_line() + facet_wrap(~property, scales = "free")
```

And for transformed data and locations
```{r}
ggplot(resid_plot_data_TP_1, aes(sample = .resid)) + stat_qq() + stat_qq_line() + facet_wrap(~property, scales = "free")
```


## Look at simulated values of distributions of random effects.
First for models fit to raw data. Note, this code plots ± 1SD, whereas merTools::plotREsim plots 95% confidence intervals by default.
```{r}
data.sc_r <- data.sc_r %>%
  mutate(randoms = map(mm_vsris, merTools::REsim))
randoms_unnest <- unnest(data.sc_r, randoms)

ggplot(filter(randoms_unnest, term == "(Intercept)"), aes(x = groupID, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd)) +
  geom_hline(yintercept=0, colour = "red") +
  facet_wrap(~property, scales = "free_y")

ggplot(filter(randoms_unnest, term == "dvlocmm"), aes(x = groupID, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd)) +
  geom_hline(yintercept=0, colour = "red") +
  facet_wrap(~property, scales = "free_y")
```

Next for models fit to transformed measurements and raw dvlocmm.
```{r}
data.sc_TP_r <- data.sc_TP_r %>%
  mutate(randoms = map(mm_vsris, merTools::REsim))
randoms_TP_unnest <- unnest(data.sc_TP_r, randoms)

ggplot(filter(randoms_TP_unnest, term == "(Intercept)"), aes(x = groupID, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd)) +
  geom_hline(yintercept=0, colour = "red") +
  facet_wrap(~property, scales = "free_y")

ggplot(filter(randoms_TP_unnest, term == "dvlocmm"), aes(x = groupID, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd)) +
  geom_hline(yintercept=0, colour = "red") +
  facet_wrap(~property, scales = "free_y")
```

Now plot for models with measured properties and dvlocmm transformed.
```{r}
data.sc_TP_r <- data.sc_TP_r %>%
  mutate(randoms_1 = map(mm_vsris_1, merTools::REsim))
randoms_TP_unnest_1 <- unnest(data.sc_TP_r, randoms_1)

ggplot(filter(randoms_TP_unnest_1, term == "(Intercept)"), aes(x = groupID, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd)) +
  geom_hline(yintercept=0, colour = "red") +
  facet_wrap(~property, scales = "free_y")

ggplot(filter(randoms_TP_unnest_1, term == "dvlocmm1"), aes(x = groupID, y = mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd)) +
  geom_hline(yintercept=0, colour = "red") +
  facet_wrap(~property, scales = "free_y")
```

Tidy up enviornment by removing unneeded data frames.
```{r}
rm(randoms_unnest, randoms_TP_unnest, randoms_TP_unnest_1, resid_plot_data, resid_plot_data_TP, resid_plot_data_TP_1)
```