1. Loading libraries
# # Uncomment below and run if these packages haven't already been installed in your version of R
# install.packages("ggplot2")
# install.packages("interactions)
# install.packages("jtools")
# install.packages("rstatix")
# install.packages("ggpubr")
# install.packages("lme4")
# install.packages("lmtest")
# install.packages("emmeans")
# load libraries
library(ggplot2) # to use ggplot for plotting
library(interactions) # to make interaction plots
library(jtools) # to use summ for summarizing models
library(rstatix) # for anova_test and get_anova_table
library(ggpubr) # for visualization, including checking qq-plots for
library(lme4) # to use lmer for mixed models (testing between and within effects)
library(lmtest) # to use bptest for testing heteroscedasticity
library(lmerTest) # to use anova for checking results
library(emmeans) # to get estimated marginal means for models
# Setting working directory path [insert your own path here]
setwd("/Users/emurbanw/OneDrive-MichiganMedicine/OneDriveDocuments/EFDC/Consulting/Jay_Kayser")
# Opening the data I created in the "Creating simulated data" chunk; making cohort a "factor" variable
data <- read.csv("sim_data.csv", header = TRUE)
data$cohort <- as.factor(data$cohort)
long <- read.csv("sim_data_long.csv", header = TRUE)
long$cohort <- as.factor(long$cohort)
2. Creating simulated data
# # I've commented this out since we only need to do this one time
#
# # Creating two dataframes for separate cohorts (cohort 1 has N = 75, cohort 2 has N = 25)
# data1 <- data.frame('ID' = 1:100, 'cohort' = 1, 'dep.1' = rnorm(n=100, mean=5.25, sd=2), 'dep.2' = rnorm(n=100, mean=3.25, sd=1), 'age' = round(rnorm(n=100, mean = 65, sd = 10)))
# data2 <- data.frame('ID' = 101:125, 'cohort' = 2, 'dep.1' = rnorm(n=25, mean=5, sd=1), 'dep.2' = rnorm(n=25, mean=2.5, sd=1), 'age' = round(rnorm(n=25, mean = 65, sd = 10)))
#
# # Combining the two dataframes
# data <- rbind(data1, data2)
#
# # Identifying factor variables
# data$cohort <- as.factor(data$cohort)
#
# # Reshaping the data into long format (necessary for lmer)
# long <- reshape(data=data,
# varying = c("dep.1", "dep.2"),
# timevar=c("time"),
# idvar=c("ID"),
# direction="long", sep=".", drop="na")
#
# # Re-ordering data by ID number
# long <- long[order(long$ID), ]
#
# # Saving the wide simulated data
# write.csv(data, "sim_data.csv", row.names = FALSE)
#
# # Saving the long simulated data
# write.csv(long, "sim_data_long.csv", row.names = FALSE)
3. Visualizing data
aggregate(dep.1 ~ cohort, data = data, mean)
## cohort dep.1
## 1 1 5.290355
## 2 2 5.225078
aggregate(dep.2 ~ cohort, data = data, mean)
## cohort dep.2
## 1 1 3.282309
## 2 2 2.333224
# This gives a plot of the mean depression score for all participants with separate lines for cohort
with(long, {interaction.plot(time, cohort, dep,
ylim = c(-1 , 10), lty = c(1, 12), lwd = 3,
ylab = "mean of dep", xlab = "time", trace.label = "cohort")
axis(side = 2, at=seq(from = -1, to = 10, by = 1), labels=seq(from = -1, to = 10, by = 1))
grid(nx = NULL, ny = 12)
})
# This gives a plot of each participant's data for the selected variables
# To view ALL participants, get rid of the [which(long$ID <=10),] bit
ggplot(long[which(long$ID <=10),]) +
aes(x = time, y = dep) +
# scale_x_discrete(breaks=c(1, 2), labels = c("pre", "post")) +
#scale_x_continuous(breaks=seq(1,2,by=1)) +
scale_x_continuous(limits=c(0.8, 2.2),breaks=c(1, 2), labels = c("pre", "post")) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
xlab("Timepoint") +
ylab("Depression score") +
facet_wrap(~ID)
4. Between-cohort comparisons
# Comparing depression ratings between the two cohorts on one timepoint
# First testing to see if variances are equal between the two cohorts at each timepoint; this determines if var.equal should be set to TRUE or FALSE
# var.equal = TRUE tells us that our test for equality of variances was supported (i.e., the test was non-significant)
# If the var.test was significant, change this to var.equal = FALSE
# paired = FALSE tells us these are independent samples
# Time point 1
var.test(data$dep.1 ~ data$cohort) # The p-value is < .05, so the variances are equal; set var.equal = FALSE in next test
##
## F test to compare two variances
##
## data: data$dep.1 by data$cohort
## F = 3.7331, num df = 99, denom df = 24, p-value = 0.0005111
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 1.843398 6.665013
## sample estimates:
## ratio of variances
## 3.733134
t.test(data$dep.1 ~ data$cohort, var.equal = FALSE, paired = FALSE)
##
## Welch Two Sample t-test
##
## data: data$dep.1 by data$cohort
## t = 0.24146, df = 74.063, p-value = 0.8099
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -0.4733924 0.6039480
## sample estimates:
## mean in group 1 mean in group 2
## 5.290355 5.225078
# Time point 2
var.test(data$dep.2 ~ data$cohort) # The p-value is > .05, so the variances are equal; set var.equal = TRUE in next test
##
## F test to compare two variances
##
## data: data$dep.2 by data$cohort
## F = 0.64154, num df = 99, denom df = 24, p-value = 0.134
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3167893 1.1453873
## sample estimates:
## ratio of variances
## 0.6415418
t.test(data$dep.2 ~ data$cohort, var.equal = TRUE, paired = FALSE)
##
## Two Sample t-test
##
## data: data$dep.2 by data$cohort
## t = 4.7854, df = 123, p-value = 4.797e-06
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## 0.5565055 1.3416639
## sample estimates:
## mean in group 1 mean in group 2
## 3.282309 2.333224
5. Within-cohort comparisons
# Assessing change over time w/ t.test (within subjects; can be done combining cohorts or run separately for each cohort; I have done both below)
# This doesn't allow a comparison of whether change over time varies as a function of being in cohort 1 or cohort 2
# Because this is paired, equality of variances is not necessary
# paired = TRUE tells us this is a paired samples t-test
# Cohort 1
t.test(data$dep.1[data$cohort==1], data$dep.2[data$cohort==1], paired = TRUE)
##
## Paired t-test
##
## data: data$dep.1[data$cohort == 1] and data$dep.2[data$cohort == 1]
## t = 10.093, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 1.613296 2.402797
## sample estimates:
## mean difference
## 2.008047
# Cohort 2
t.test(data$dep.1[data$cohort==2], data$dep.2[data$cohort==2], paired = TRUE)
##
## Paired t-test
##
## data: data$dep.1[data$cohort == 2] and data$dep.2[data$cohort == 2]
## t = 9.4472, df = 24, p-value = 1.475e-09
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 2.260077 3.523630
## sample estimates:
## mean difference
## 2.891854
# Both cohorts together
t.test(data$dep.1, data$dep.2, paired = TRUE)
##
## Paired t-test
##
## data: data$dep.1 and data$dep.2
## t = 12.632, df = 124, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 1.842475 2.527142
## sample estimates:
## mean difference
## 2.184808
6. Checking ANOVA assumptions
# Identifying outliers
long %>%
group_by(time, cohort) %>%
identify_outliers(dep)
## # A tibble: 4 × 7
## cohort time ID age dep is.outlier is.extreme
## <fct> <int> <int> <int> <dbl> <lgl> <lgl>
## 1 1 1 6 62 10.7 TRUE FALSE
## 2 1 1 54 86 9.92 TRUE FALSE
## 3 1 1 86 65 -0.554 TRUE FALSE
## 4 1 2 32 71 5.60 TRUE FALSE
# Checking normality of distribution; p-values should be > .05
# Shapiro test is sensitive to sample sizes over 50, so even small deviations may lead to significant results
long %>%
group_by(time, cohort) %>%
shapiro_test(dep)
## # A tibble: 4 × 5
## cohort time variable statistic p
## <fct> <int> <chr> <dbl> <dbl>
## 1 1 1 dep 0.990 0.636
## 2 2 1 dep 0.993 0.999
## 3 1 2 dep 0.996 0.994
## 4 2 2 dep 0.949 0.234
# Checking homogeneity of variances; p-values should be > .05
long %>%
group_by(time) %>%
levene_test(dep ~ cohort)
## # A tibble: 2 × 5
## time df1 df2 statistic p
## <int> <int> <int> <dbl> <dbl>
## 1 1 1 123 6.88 0.00982
## 2 2 1 123 1.56 0.214
# Checking assumption of sphericity
# This is tested and corrected for automatically when running anova_test and get_anova_table
# Checking homogeneity of covariances; p-values should be > .001; sensitive to unequal sample sizes
box_m(long[, "dep", drop = FALSE], long$cohort)
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 0.000548 0.981 1 Box's M-test for Homogeneity of Covariance Matric…
7. Between/Within cohort comparisons
# THERE ARE SEVERAL OPTIONS HERE, AND IT DEPENDS ON CHARACTERISTICS OF YOUR DESIGN
# For example, if you have unbalanced cohort sizes or missing data, you probably will need to run the lmer command which is more flexible
# Otherwise, if you have a very simple design you could probably stick with a repeated measures or mixed GLM.
# Often they will give you similar results in terms of rejecting or failing to reject the null hypothesis, but not always.
# In the example here, there are actually discrepancies in the significance of the cohort term, even though we don't necessarily want to interpret that anyway in the presence of an interaction.
# So for this example dataset, I recommend using model2 given the unbalanced design, presence of outliers, and flexibility of lmer.
########## Effect of time in a Repeated Measures-GLM ##########
model0 <- anova_test(data = long, dv = dep, wid = ID, within = time)
# Checking the results
get_anova_table(model0)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 time 1 124 159.567 5.06e-24 * 0.38
########## Effect of time and cohort in a Mixed-GLM ##########
model1 <- anova_test(data = long, dv = dep, wid = ID, within = time, between = cohort)
# Checking the results
get_anova_table(model1)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 cohort 1 123 5.172 2.50e-02 * 0.021
## 2 time 1 123 131.819 3.48e-21 * 0.339
## 3 cohort:time 1 123 4.289 4.00e-02 * 0.016
########## Linear mixed effects model w/ random intercept ##########
# Running the model; (1|ID) specifies that observations are nested within ID and that intercepts of the model can vary between individuals
model2 <- lmer(dep ~ time*cohort + (1|ID), data = long)
# Checking the results
anova(model2)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## time 240.090 240.090 1 123.00 131.8191 < 2e-16 ***
## cohort 2.655 2.655 1 152.42 1.4579 0.22913
## time:cohort 7.811 7.811 1 123.00 4.2886 0.04046 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summ(model2, digits=3, confint=TRUE, part.corr = TRUE, robust=FALSE)
## MODEL INFO:
## Observations: 250
## Dependent Variable: dep
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 884.128, BIC = 905.257
## Pseudo-R² (fixed effects) = 0.400
## Pseudo-R² (total) = 0.427
##
## FIXED EFFECTS:
## -------------------------------------------------------------------------
## Est. 2.5% 97.5% t val. d.f. p
## ------------------ -------- -------- -------- --------- --------- -------
## (Intercept) 7.298 6.704 7.893 24.074 152.424 0.000
## time -2.008 -2.382 -1.634 -10.521 123.000 0.000
## cohort2 0.819 -0.510 2.147 1.207 152.424 0.229
## time:cohort2 -0.884 -1.720 -0.047 -2.071 123.000 0.040
## -------------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ------------------------------------
## Group Parameter Std. Dev.
## ---------- ------------- -----------
## ID (Intercept) 0.290
## Residual 1.350
## ------------------------------------
##
## Grouping variables:
## --------------------------
## Group # groups ICC
## ------- ---------- -------
## ID 125 0.044
## --------------------------
summary(model2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dep ~ time * cohort + (1 | ID)
## Data: long
##
## REML criterion at convergence: 872.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0975 -0.5115 -0.0109 0.4929 3.8494
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.0841 0.29
## Residual 1.8214 1.35
## Number of obs: 250, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.2984 0.3032 152.4240 24.074 <2e-16 ***
## time -2.0080 0.1909 123.0000 -10.521 <2e-16 ***
## cohort2 0.8185 0.6779 152.4240 1.207 0.2291
## time:cohort2 -0.8838 0.4268 123.0000 -2.071 0.0405 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time cohrt2
## time -0.944
## cohort2 -0.447 0.422
## time:cohrt2 0.422 -0.447 -0.944
# Looking at estimated marginal means based on the model
model2_means <- emmeans(model2, ~time*cohort)
summary(model2_means)
## time cohort emmean SE df lower.CL upper.CL
## 1 1 5.29 0.138 246 5.02 5.56
## 2 1 3.28 0.138 246 3.01 3.55
## 1 2 5.23 0.276 246 4.68 5.77
## 2 2 2.33 0.276 246 1.79 2.88
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
interact_plot(model2, pred = time, modx = cohort, plot.points = TRUE, partial.residuals= FALSE, interval = FALSE, facet.modx = FALSE, pred.labels = c("pre", "post"), x.label = "Timepoint", y.label = "Depression Score", legend.main = "Cohort")
########## Linear mixed effects model w/ random intercept and slope ##########
# # The below model would add that change in depression scores over time can vary between individuals
# # However, this will not work with only two timepoints, as it overparamaterizes the model by trying to estimate random intercepts AND random slopes.
# model3 <- lmer(dep ~ time*cohort + (1 + time|ID), data = long)
# anova(model3)
# summ(model3, digits=3, confint=TRUE, part.corr = TRUE, robust=FALSE)
# summary(model3)