# ダウンロードデータの説明 # 診断と治療社『保健医療従事者のためのマルチレベル分析活用ナビ』より # http://www.shindan.co.jp/download/index.php?pcode=205300 x <- read.csv("http://www.shindan.co.jp/download/205300/cholesterol.csv") # graph1 layout(matrix(1:4, 2, byrow=TRUE)) hist(x$cholesterol) hist(x$cholesterol_base) matplot(rbind(x$cholesterol_base, x$cholesterol), type="l", col=x$intervention+2, lty=1, lwd=1, axes=FALSE) axis(1, 1:2, c("base","after")) axis(2, 3:7*50) stripchart(cholesterol-cholesterol_base ~ intervention, data=x, method="stack", ylab="intervention") # t-test t.test(x$cholesterol_base, x$cholesterol, paired=TRUE, var.equal=FALSE) # graph2 layout(1) plot(cholesterol ~ cholesterol_base, pch=intervention+1, data=x) legend("topleft", pch=1:2, legend=c("without intervention","with intervention")) # ANCOVA x$interventionF <- as.factor(x$intervention) res1 <- lm(cholesterol ~ cholesterol_base * interventionF, data=x) summary(res1) res2 <- lm(cholesterol ~ cholesterol_base + interventionF, data=x) summary(res2) # regressions separately # (natural for this data, but unsuitable for the purpose) summary(lm(cholesterol ~ cholesterol_base, data=subset(x, intervention==0))) summary(lm(cholesterol ~ cholesterol_base, data=subset(x, intervention==1))) # graph3 plot(cholesterol ~ cholesterol_base, pch=intervention+1, data=x, col=topo.colors(10)[clinic]) legend("topleft", pch=1:2, legend=c("without intervention","with intervention")) legend("bottomright", pch=2, col=topo.colors(10), legend=1:10, title="clinic") # multilevel library(lmerTest) res3 <- lmer(cholesterol ~ cholesterol_base + intervention + (1 | clinic), data=x, REML=FALSE) summary(res3) confint(res3) # res4 <- lm(cholesterol ~ cholesterol_base + intervention, # data=x) # summary(res4) # anova(res4, res3) might make likelihood ratio test? # But it's impossible. And, lmer cannot apply for # the model without random effect.