data set can be seen through this link: https://uoepsy.github.io/data/dapr3_2122_report1.csv.
The data set has information on participants from different hospitals who have had a stroke and is trying to asses whether their birthweight affects the participants change in cognitive function following the stoke. Cognitive function is assessed (straight after the stroke, 6,12,18 and 24 months following the stoke).
This is what iv done so far
library(tidyverse)
library(car)
library(Matrix)
library(lme4)
# Read in data set
data <- read_csv("https://uoepsy.github.io/data/dapr3_2122_report1.csv")
head(data)
# Remove N/A from data set
newdata <- na.omit(data)
#LM of birthweight effect on ACER when control age + educ
lm <- lm(ACER ~ age + educ, data = newdata)
summary(lm)
# Cooks Distance
cooksd = cooks.distance(lm)
influential <- as.numeric(names(cooksd) [(cooksd > 4/694)])
cleandata <- newdata[-influential,]
# Null model
m.null <- lmer(ACER ~ 1 + (1 | pid), data=cleandata)
summary(m.null)
We can see the 0.35 / (0.35 + 6.54), of the total variance is attributable to participant-level variation.
#Now lets suppose we want to compare this null model with a model with an effect of birthweight(to assess whether there is overall change over time).
modA <- lmer(ACER ~ 1 + birthweight + (1 + birthweight | pid), data=cleandata)
modB <- lmer(ACER ~ 1 + birthweight + (1 | pid), data=cleandata)
isSingular(modA)
isSingular(modB)
#If we want to conduct a model comparison to isolate the effect of overall change over time (a fixed effect of birthweight), we might want to compare these two models:
m.base0 <- lmer(ACER ~ 1 + (1 + birthweight | pid), data=cleandata)
m.base <- lmer(ACER ~ 1 + birthweight + (1 + birthweight | pid), data=cleandata)
isSingular(m.base0)
isSingular(m.base)
# Straightforward LRT
anova(m.base0, m.base)
# parametric bootstrap LRT
library(pbkrtest)
PBmodcomp(m.base, m.base0)