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 controlling for age and education
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,]
cleandata %>%
mutate(hospppt = interaction(hospital,pid)) %>%
summarise(across(everything(), n_distinct))
L3_mdl <- cleandata %>% mutate(across(c(birthweight, hospital, pid, age, educ)))
L3_mdl0 <- lmer(ACER ~ timepoint * birthweight + (1 | hospital/pid),data = cleandata)
plot(L3_mdl0, type=c("p","smooth"))
Compared to the model (black lines bellow), the participents scores (coloured lines) are more
library(broom.mixed)
augment(L3_mdl0) %>%
ggplot(aes(x=timepoint, col=pid)) +
geom_point(aes(y = ACER))+
geom_path(aes(y = ACER))+
geom_path(aes(y = .fitted), col="black", alpha=.3)+
guides(col="none")+
facet_wrap(~pid)
L3_mdl1 <- lmer(ACER ~ timepoint * birthweight + (1 + timepoint | hospital/pid), data = cleandata)
plot(L3_mdl1, type=c("p","smooth"))
library(HLMdiag)
infl <- hlm_influence(L3_mdl1, level=1)
dotplot_diag(infl$cooksd, cutoff = "internal")
infpid <- hlm_influence(L3_mdl1, level= "pid:hospital")
dotplot_diag(infpid$cooksd, cutoff = "internal", index = infpid$`pid:hospital`) +
scale_y_continuous(limits = c(0,.05))
`
del <- case_delete(L3_mdl1, level = "pid:hospital", delete = "10:HSP12")
cbind(del$fixef.original, del$fixef.delete)
del1 <- case_delete(L3_mdl1, level = "pid:hospital", delete = "4:HSP10")
cbind(del1$fixef.original, del1$fixef.delete)
del2 <- case_delete(L3_mdl1, level = "pid:hospital", delete = "3:HSP5")
cbind(del2$fixef.original, del2$fixef.delete)
del3 <- case_delete(L3_mdl1, level = "pid:hospital", delete = "4:HSP11")
cbind(del3$fixef.original, del3$fixef.delete)
infhosp <- hlm_influence(L3_mdl1, level = "hospital")
dotplot_diag(infhosp$cooksd, cutoff = "internal", index = infhosp$hospital)
L3_null <- lmer(ACER ~ 1 + (1 | hospital/pid), data = newdata, control = lmerControl(optimizer = "bobyqa"))
as.data.frame(VarCorr(L3_null)) %>%
select(grp, vcov) %>%
mutate(
icc = cumsum(vcov)/sum(vcov)*100
)
library(lmeresampler)
L3_mdlBS <- bootstrap(L3_mdl1, .f=fixef, type = "case", B = 2000, resample = c(FALSE, FALSE, TRUE))
is.atomic(L3_mdlBS)
summary(L3_mdlBS)
confint(L3_mdlBS, type = "perc")
That is my code above, I think I might need to dress the warnings on rescaling variables first and then the bootstrap model might work but I don't know how to rescale my variables?