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 res1 <- lm(cholesterol ~ cholesterol_base * as.factor(intervention), data=x) summary(res1) res2 <- lm(cholesterol ~ cholesterol_base + as.factor(intervention), data=x) summary(res2) # 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)