协方差分析

rm(list=ls() )
setwd("C:/Users/DELL/Desktop")
ss1 <- 37#对照组
ss2 <- 40#实验组
set.seed(203)
age1 <- round(runif(ss1,3,42),0);stats::shapiro.test(age1);mean(age1);sd(age1)
set.seed(209)
age2 <- round(runif(ss2,3,42),0);stats::shapiro.test(age2);mean(age2);sd(age2)
set.seed(203)
sex1 <- sample(c(1,2),ss1,replace = T);table(sex1)
set.seed(204)
sex2 <- sample(c(1,2),ss2,replace = T);table(sex2)
set.seed(203)
bZtr1 <- round(rnorm(ss1,65,8));stats::shapiro.test(bZtr1)
mean(bZtr1);sd(bZtr1)
set.seed(204)
aZtr1 <- bZtr1+sample(15:23,ss1,replace = T);stats::shapiro.test(aZtr1)
mean(aZtr1);sd(aZtr1)
set.seed(203)
bZtr2 <- round(rnorm(ss2,65,8));mean(bZtr2);stats::shapiro.test(bZtr2)
mean(bZtr2);sd(bZtr2)
set.seed(20201)
aZtr2 <- bZtr2+sample(18:26,ss2,replace = T);stats::shapiro.test(bZtr2)
mean(aZtr2);sd(aZtr2)

set.seed(200)
bYtr1 <- round(rnorm(ss1,65,8));stats::shapiro.test(bYtr1)
mean(bYtr1);sd(bYtr1)
set.seed(206)
aYtr1 <- bYtr1+sample(15:23,ss1,replace = T);stats::shapiro.test(aYtr1)
mean(aYtr1);sd(aYtr1)
set.seed(205)
bYtr2 <- round(rnorm(ss2,65,8));stats::shapiro.test(bYtr2)
mean(bYtr2);sd(bYtr2)
set.seed(202)
aYtr2 <- bYtr2+sample(18:26,ss2,replace = T);stats::shapiro.test(bYtr2)
mean(aYtr2);sd(aYtr2)

c    <- cbind.data.frame(age1,sex1,bZtr1,aZtr1,bYtr1,aYtr1)
c$g  <- 1
c$aZtr1[c$aZtr1>=100] <-95
c$aYtr1[c$aYtr1>=100] <-95
c    <-dplyr::rename(c,age=age1,sex=sex1,bZtr=bZtr1,aZtr=aZtr1,bYtr=bYtr1,aYtr=aYtr1)
t    <- cbind.data.frame(age2,sex2,bZtr2,aZtr2,bYtr2,aYtr2)
t$g  <- 2
t$aZtr2[t$aZtr2>=100] <-95
t$aYtr2[t$aYtr2>=100] <-95
t    <-dplyr::rename(t,age=age2,sex=sex2,bZtr=bZtr2,aZtr=aZtr2,bYtr=bYtr2,aYtr=aYtr2)
data <- rbind.data.frame(c,t)
data$g<- as.factor(data$g)
head(data)
write.csv(data,file = "data.csv",row.names = F)
#年龄
t.test(age1,age2,paired = FALSE,var.equal = T) 
#性别
sex <- table(data$g,data$age)
chisq.test(sex)
#正态性可视化
library(car)
qqPlot(lm(aZtr ~ bZtr + g , data = data), simulate = T, 
       main = "QQ Plot", labels = FALSE)

qqPlot(lm(aYtr ~ bYtr + g , data = data), simulate = T, 
       main = "QQ Plot", labels = FALSE)
#方差齐
library(stats) 
bartlett.test(aZtr ~ g, data = data)
bartlett.test(aYtr ~ g, data = data)
#斜率相等==交互无意义
library(multcomp)
fitZ <- aov(aZtr ~ bZtr * g,data = data)
summary(fitZ)
fitY <- aov(aYtr ~ bYtr * g,data = data)
summary(fitY)

#协方差
library(HH)
fitZ1 <-ancova(aZtr ~ bZtr + g, data = data);fitZ1
fitY1 <-ancova(aYtr ~ bYtr + g, data = data);fitY1
library(effects)
lsmeanZ <- effect(term='g',mod = fitZ1,se=list(compute = T, level = 0.95))
effect(term='g',mod = fitZ1,se=list(compute = T, level = 0.95)) #调整后均值
ciupperZ <- c(lsmeanZ$lower[1],lsmeanZ$upper[1]);ciupperZ
cilowerZ <- c(lsmeanZ$lower[2],lsmeanZ$upper[2]);cilowerZ

lsmeanY <- effect(term='g',mod = fitY1,se=list(compute = T, level = 0.95))
effect(term='g',mod = fitY1,se=list(compute = T, level = 0.95)) #调整后均值
ciupperY <- c(lsmeanY$lower[1],lsmeanY$upper[1]);ciupperY
cilowerY <- c(lsmeanY$lower[2],lsmeanY$upper[2]);cilowerY

library(multcomp)       
contrast <-rbind('1vs2'=c(-1,1))
diffmeanZ <-glht(fitZ1,linfct=mcp(g=contrast));diffmeanZ #调整后均值之差
diffmeanY <-glht(fitY1,linfct=mcp(g=contrast));diffmeanY #调整后均值之差

参考另一篇，结合着看