zoukankan      html  css  js  c++  java

  • R语言学习笔记（七）：方差分析

    单因素方程分析
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
with 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 p
19 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))

    用回归来做ANOVA

library(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
  • 相关阅读:
    MySQL-sql语句
    常见IP端口
    Java开发异常
    微信开发-charles抓包
    6、分组查询
    5、关联查询
    4、查询
    3、表的管理
    2、用户管理
    1、oracle基本使用
  • 原文地址：https://www.cnblogs.com/GhostBear/p/7724667.html
Copyright © 2011-2022 走看看