单因素方程分析install.packages("multcomp")library(multcomp)attach(cholesterol)table(trt)aggregate(response, by=list(trt),FUN=mean)aggregate(response, by=list(trt),FUN=sd)fit<-aov(response~trt)summary(fit) 　　     　　Df 　　Sum Sq 　　Mean Sq 　　F value 　　Pr(>F) trt 　　        4         1351.4           337.8         32.43          9.82e-13 ***Residuals    45        468.8            10.4 ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

install.packages("gplots")library(gplots)plotmeans(response~trt,xlab="Treatment",ylab="Response", main="Mean Plot\nwith 95% CI")detach(cholesterol)

#多重比较TukeyHSD(fit)par(las=2)par(mar=c(5,8,4,2))plot(TukeyHSD(fit))

 

library(multcomp)par(mar=c(5,4,6,2))tuk<-glht(fit,linfct=mcp(trt="Tukey"))plot(cld(tuk,level=.05),col="lightgrey")

 

#离群点检测 - #评估检验的假设条件library(car)qqPlot(lm(response~trt,data=cholesterol),simulate=TRUE,main="Q-Q Plot",labels=FALSE)

#Bartlett检验bartlett.test(response~trt,data=cholesterol)

Bartlett test of homogeneity of variances

data: response by trt

Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653

outlierTest(fit) #离群点检测No Studentized residuals with Bonferonni p < 0.05  没发现离群点Largest |rstudent|:rstudent unadjusted p-value Bonferonni p19 2.251149 0.029422 NA

 

单因素协方差分析data(litter,package="multcomp")attach(litter)table(dose)weight<-weight[c(1:60)]gesttime<-gesttime[c(1:60)]aggregate(weight,by=list(dose),FUN=mean)fit<-aov(weight~gesttime+dose)summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

gesttime 1 107.0 107.04 7.099 0.0101 *

dose 2 45.9 22.97 1.523 0.2269

Residuals 56 844.3 15.08

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(effects)effect("dose",fit)

dose effect

dose

0.5 1 2

30.44530 31.68956 29.46164

 

#多重比较library(multcomp)contrast<-rbind("no drug vs. drug"=c(3,-1,-1,-1))summary(glht(fit,linfct=mcp(dose=contrast)))#评估验证的假设条件library(multcomp)fit2<-aov(weight~gesttime*dose,data=litter)summary(fit2)#结果可视化install.packages("HH")library(HH)ancova(weight~gesttime+dose,data=litter)

 

双因素分析attach(ToothGrowth)table(supp,dose)aggregate(len,by=list(supp,dose),FUN=mean)aggregate(len,by=list(supp,dose),FUN=sd)dose<-factor(dose)fit<-aov(len~supp*dose)summary(fit)

Df Sum Sq Mean Sq F value Pr(>F)

supp 1 205.4 205.4 15.572 0.000231 ***

dose 2 2426.4 1213.2 92.000 < 2e-16 ***

supp:dose 2 108.3 54.2 4.107 0.021860 *

Residuals 54 712.1 13.2

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

 

interaction.plot(dose,supp,len,type="b",col=c("red","blue"),pch=c(16,18),main="Interaction between Dose and Supplement Type")

detach(ToothGrowth)

 

library(gplots)plotmeans(len~interaction(supp,dose,sep=" "),connect = list(c(1,3,5),c(2,4,6)),col=c("red","darkgreen"),main="Interaction Plot with 95% CIs",xlab="Treatment and Dose Combination")

library(HH)interaction2wt(len~supp*dose)

重复测量方差分析CO2$conc<-factor(CO2$conc)w1b1<-subset(CO2,Treatment=='chilled')fit<-aov(uptake~conc*Type+Error(Plant/(conc)),w1b1)summary(fit)

Error: Plant

Df Sum Sq Mean Sq F value Pr(>F)

Type 1 2667.2 2667.2 60.41 0.00148 **

Residuals 4 176.6 44.1

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Plant:conc

Df Sum Sq Mean Sq F value Pr(>F)

conc 6 1472.4 245.40 52.52 1.26e-12 ***

conc:Type 6 428.8 71.47 15.30 3.75e-07 ***

Residuals 24 112.1 4.67

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

 

par(las=2)

par(mar=c(10,4,4,2))

with(w1b1,interaction.plot(conc,Type,uptake,type="b",col=c("red","blue"),pch=c(16,18),main="Interaction Plot for Plant Type and Concentration"))

boxplot(uptake~Type*conc, data=w1b1,col=(c("gold","green")),main="Chilled Quebec and Mississippi Plants",ylab="Carbon dioxide uptake rate (umol/m^2 sec)")

 

多元方差分析library(MASS)attach(UScereal)shelf<-factor(shelf)y<-cbind(calories,fat,sugars)aggregate(y,by=list(shelf),FUN=mean)cov(y)

calories fat sugars

calories 3895.24210 60.674383 180.380317

fat 60.67438 2.713399 3.995474

sugars 180.38032 3.995474 34.050018

fit<-manova(y~shelf) summary(fit)

Df Pillai approx F num Df den Df Pr(>F)

shelf 2 0.4021 5.1167 6 122 0.0001015 ***

Residuals 62

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary.aov(fit)

Response calories :

Df Sum Sq Mean Sq F value Pr(>F)

shelf 2 50435 25217.6 7.8623 0.0009054 ***

Residuals 62 198860 3207.4

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Response fat :

Df Sum Sq Mean Sq F value Pr(>F)

shelf 2 18.44 9.2199 3.6828 0.03081 *

Residuals 62 155.22 2.5035

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Response sugars :

Df Sum Sq Mean Sq F value Pr(>F)

shelf 2 381.33 190.667 6.5752 0.002572 **

Residuals 62 1797.87 28.998

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

评估假设检验center<-colMeans(y)n<-nrow(y)p<-ncol(y)cov<-cov(y)d<-mahalanobis(y,center,cov)coord<-qqplot(qchisq(ppoints(n),df=p),d,main="Q-Q Plot Assessing Multivariate Normality",ylab="Mahalanobis D2")abline(a=0,b=1)identify(coord$x,coord$y,labels=row.names(UScereal))

 

 

用回归来做ANOVAlibrary(multcomp)levels(cholesterol$trt)[1] "1time"  "2times" "4times" "drugD"  "drugE"  fit.aov<-aov(response~trt,data=cholesterol)summary(fit.aov)

Df Sum Sq Mean Sq F value Pr(>F)

trt 4 1351.4 337.8 32.43 9.82e-13 ***

Residuals 45 468.8 10.4

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.lm<-lm(response~trt,data=cholesterol)summary(fit.lm)

Call:

lm(formula = response ~ trt, data = cholesterol)

Residuals:

Min 1Q Median 3Q Max

-6.5418 -1.9672 -0.0016 1.8901 6.6008

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 5.782 1.021 5.665 9.78e-07 ***

trt2times 3.443 1.443 2.385 0.0213 *

trt4times 6.593 1.443 4.568 3.82e-05 ***

trtdrugD 9.579 1.443 6.637 3.53e-08 ***

trtdrugE 15.166 1.443 10.507 1.08e-13 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.227 on 45 degrees of freedom

Multiple R-squared: 0.7425, Adjusted R-squared: 0.7196

F-statistic: 32.43 on 4 and 45 DF, p-value: 9.819e-13

contrasts(cholesterol$trt)

2times 4times drugD drugE

1time 0 0 0 0

2times 1 0 0 0

4times 0 1 0 0

drugD 0 0 1 0

drugE 0 0 0 1

fit.lm<-lm(response~trt,data=cholesterol,contrasts="contr.helmert")summary(fit.lm)

Call:

lm(formula = response ~ trt, data = cholesterol, contrasts = "contr.helmert")

Residuals:

Min 1Q Median 3Q Max

-6.5418 -1.9672 -0.0016 1.8901 6.6008

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 5.782 1.021 5.665 9.78e-07 ***

trt2times 3.443 1.443 2.385 0.0213 *

trt4times 6.593 1.443 4.568 3.82e-05 ***

trtdrugD 9.579 1.443 6.637 3.53e-08 ***

trtdrugE 15.166 1.443 10.507 1.08e-13 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.227 on 45 degrees of freedom

Multiple R-squared: 0.7425, Adjusted R-squared: 0.7196

F-statistic: 32.43 on 4 and 45 DF, p-value: 9.819e-13