生物统计

《白话统计》读书笔记-t and anovar

2019-04-19  本文已影响39人  drlee_fc74

t检验

t检验的使用条件

  1. t检验用于两组比较。如果是多组比较则不使用
  2. 如果数据严重偏态,建议不要使用t检验。但是如果数据稍微有些偏态那么也是可以使用的。如果偏态则可以考虑:变量转换或者可是使用Wilcoxon秩和
  3. 如何两组方差相差较大,则不建议使用t检验。如果方差大可以考虑:变量转换或者可是使用Wilcoxon秩和

t检验的实现

独立样本t检验

独立样本t检验使用的函数为t.test(y ~ x, data)其中,y是连续性的,x是二分类变量。或者可是使用t.test(y1, y2)两者都是连续性的变量。

high <- c(134, 146, 104, 119, 124, 161, 107, 83, 113, 129, 97, 123)
low <- c(70, 118, 101, 85, 107, 132, 94)
t.test(high, low, paired = F)
## 
##  Welch Two Sample t-test
## 
## data:  high and low
## t = 1.9107, df = 13.082, p-value = 0.07821
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.469073 40.469073
## sample estimates:
## mean of x mean of y 
##       120       101

配对样本t检验

相较于独立样本,配对样本t检验只是把paired更为T即可

#输入两组值
x <- c(2.41,2.90,2.75,2.23,3.67,4.49,5.16,5.45,2.06,1.64,1.06,0.77)
y <- c(2.80,3.04,1.88,3.43,3.81,4.00,4.44,5.41,1.24,1.83,1.45,0.92)
#配对样本t检验
t.test(x,y,paired=T)
## 
##  Paired t-test
## 
## data:  x and y
## t = 0.16232, df = 11, p-value = 0.874
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3558501  0.4125167
## sample estimates:
## mean of the differences 
##              0.02833333

Wilcoxon秩和的实现

Wilcoxon秩和检验可以通过wilcox.test来实现.同样的独立和配对样本之间的区别只在于pairedT/F

wilcox.test(high, low)
## Warning in wilcox.test.default(high, low): cannot compute exact p-value
## with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  high and low
## W = 62.5, p-value = 0.09083
## alternative hypothesis: true location shift is not equal to 0

方差分析

R当中,如果需要进行方差分析可以使用aov(formula,data)函数。

方差分析中变异分解的思想

我们在进行多组比较的时候,比较不同组之间的差异,应该考虑两点内容:个体之间的差异以及分组之间的差异。一般来说我们就是来比较这两个内容到底是那个主要影响变异的产生。因此我们需要做的就是:求出组间方差和组内方差比较两者是否有差异以及差异变化是什么样的即可。 组间方差/组内方差=F值.如果F很大则代表组间变异大,如果F值小则代表组内变异大。

方差分析中的F和线性回归当中的F

我们在线性回归当中也是可以F值的。其实通过上面的理解F值代表说不同分组之间的变异和其他变异的差异。同样的,线性回归也可以理解为,我们找到的模型对于整个数据的变异和其他变异的差异。因此其实方差和线性回归当中都出现F值也是可以理解的。只不过对于结局变量来说,方差是分类变量,而线性回归则是连续性变量 当然,其实回归分析当中还有类似R^2来评价变异的

fit <- lm(disp ~ hp, data = mtcars)
summary(fit)
## 
## Call:
## lm(formula = disp ~ hp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -198.97  -37.39  -18.27   57.55  157.91 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.9925    32.6066   0.644    0.525    
## hp            1.4298     0.2019   7.080 7.14e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 77.09 on 30 degrees of freedom
## Multiple R-squared:  0.6256, Adjusted R-squared:  0.6131 
## F-statistic: 50.13 on 1 and 30 DF,  p-value: 7.143e-08

模型结果的最后则是F值。

基于方差分析的不同实验设计

单因素方差分析

这类实验设计只包括一个考虑因素。例如:不同的药物对于治疗疗效的评价。

library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
data("cholesterol")
str(cholesterol)
## 'data.frame':    50 obs. of  2 variables:
##  $ trt     : Factor w/ 5 levels "1time","2times",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ response: num  3.86 10.39 5.91 3.06 7.72 ...
fit <- aov(response ~ trt, data = cholesterol)
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

析因设计的方差分析

head(ToothGrowth)
##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5
fit <- aov(len ~ supp*dose, data = ToothGrowth)
summary(fit)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## supp         1  205.4   205.4  12.317 0.000894 ***
## dose         1 2224.3  2224.3 133.415  < 2e-16 ***
## supp:dose    1   88.9    88.9   5.333 0.024631 *  
## Residuals   56  933.6    16.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果当中supp:dose代表两者交互的结果。

library(tidyverse)
## ─ Attaching packages ────────────────── tidyverse 1.2.1 ─
## ✔ ggplot2 3.1.0       ✔ purrr   0.3.2  
## ✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0
## ─ Conflicts ─────────────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
ToothGrowth %>% group_by(dose, supp) %>% summarise(mean = mean(len)) %>% 
    ggplot(., aes(factor(dose), mean, color = supp, group = supp)) + geom_point() + geom_line()
image.png
HH::interaction2wt(len~supp*dose, data = ToothGrowth)
image.png

多元方差分析

当因变量不止一个时,可以使用多元方差分析对他们进行同时分析。 我们可以使用manova来进行分析

library(MASS)
data("UScereal")
smalldat <- UScereal[,c("calories", "fat", "sugars", "shelf")]
head(smalldat)
##                           calories      fat   sugars shelf
## 100% Bran                 212.1212 3.030303 18.18182     3
## All-Bran                  212.1212 3.030303 15.15151     3
## All-Bran with Extra Fiber 100.0000 0.000000  0.00000     3
## Apple Cinnamon Cheerios   146.6667 2.666667 13.33333     1
## Apple Jacks               110.0000 0.000000 14.00000     2
## Basic 4                   173.3333 2.666667 10.66667     3
fit <- manova(cbind(calories, fat, sugars) ~ shelf, data = smalldat)
summary.aov(fit)  ###输出单变量结果
##  Response calories :
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## shelf        1  45313   45313  13.995 0.0003983 ***
## Residuals   63 203982    3238                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response fat :
##             Df  Sum Sq Mean Sq F value   Pr(>F)   
## shelf        1  18.421 18.4214   7.476 0.008108 **
## Residuals   63 155.236  2.4641                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response sugars :
##             Df  Sum Sq Mean Sq F value  Pr(>F)  
## shelf        1  183.34  183.34   5.787 0.01909 *
## Residuals   63 1995.87   31.68                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

方差分析后的两两比较方法选择

多于2组的分析是要用方差分析的。统一分组完之后,如果还要进行单个的两两比较。那么就需要选择不同的方式了。 1. 如果比较的组数较多(4组以上),如果想要进行两两比较,可以选择Turkey或者REGWQ法。 2. 如果组数不多(3组),则可以选择Turkey或者Bonferroni法 3. 如果不仅想两两比较,还行进行其他方式的比较,如一个对照组和另外两个实验组的均值比较。则建议选择BonferroniScheffe法。如果组数较多则不建议选择Bonferroni 4. 如果有一个明确的对照组和多个实验组,分别比较各实验组和对照组,可以选用Dunnett法。 5. 如果各组例数相差较大,尤其是相差明显,可以考虑Games-Howell法。 6. 如果不是进行组间比较,而是进行多指标筛选,那么可以考虑FDR/Holm/Hochberg

方差分析的替代—Kruskal-Wallis秩和检验

如果不正态或者方差不齐的话,可以通过Kruskal-Wallis检验来进行替代。 R语言中可以使用kruskal.test来进行检验

fit <- kruskal.test(fat ~ shelf, data = smalldat)
fit
## 
##  Kruskal-Wallis rank sum test
## 
## data:  fat by shelf
## Kruskal-Wallis chi-squared = 5.5524, df = 2, p-value = 0.06228

秩和当中的两两比较

秩和当中的两两比较有三种方式分别是DSCF/Conover-lman/Nemenyi法。 R语言中可以使用PMCMRplus包中的posthoc.kruskal.dunn.testposthoc.kruskal.conover.test;以及posthoc.kruskal.nemeyi.test来分别实现三种统计。一般来说如果是进行三组或者四组两两比较的话,可以使用Nemenyi/DSCF法。

上一篇下一篇

猜你喜欢

热点阅读