R常见统计模型使用

数据描述：

summary(A) psych::describe(A)

直方图：

hist(A)

概率密度图density plot:

hist(A,prob=TRUE) lines(density(A, bw="SJ"))

累计分布图CDF：

plot(ecdf(A), do.points=FALSE, verticals=TRUE)

拟合同参数高斯分布的CDF：

lines(x, pnorm(x, mean=mean(A), sd=sqrt(var(A))), lty=3)

双样本检验

经典的two-sample t-test：

t.test(A, B, var.equal=TRUE)

如果不满足同方差，则用 Welch modified two-sample t-test :

res=t.test(A, B) res$p.val

检验同方差的方法可用F-test:

var.test(A, B)

如果不满足正态性，则用非参数检验，即two-sample Mann-Whitney test:（但是仍然assume两个分布是相似的）

wilcox.test(A, B)

如果连分布也很不相似，则考虑用permutation test

coin::independence_test(Length ~ Hand, data = Data)

是否满足正态分布，有三种方法检验：

1. QQ图：

    par(pty="s") # arrange for a square figure region qqnorm(long) qqline(long)

2. Shapiro-Wilk test：

    shapiro.test(A) # Shapiro-Wilk normality test

3. Kolmogorov-Smirnov test：

    ks.test(A, "pnorm", mean = mean(A), sd = sqrt(var(A)))

独立性检验：

卡方分布有一个重要应用就是根据样本数据来检验两个分类变量的独立性，我们以CO2数据为例来说明chisq.test函数的使用，help(CO2)可以了解更多信息。

data(CO2) #读入内置的数据包，其中Type和Treatmen是其中两个分类变量。 
chisq.test(table(CO2$Type,CO2$Treatment)) #检验这两个因子之间是否独立

结果显示P值为0.82，因此可以认为两因子之间独立。在样本较小的情况下，还可以使用fisher精确检验，对应的函数是fisher.test。

proportion test:(在两个group的情况下，结果和卡方检验是一样的)

library(MASS) 
table(quine$Eth, quine$Sex) 
prop.test(table(quine$Eth, quine$Sex), correct=FALSE)

ANOVA：

anova(lm(y~x)) 
summary(aov(y~x)) # summary(aov(lm(y~x)))也是一样的结果 
summary(aov(y~x+covariate)) # ANCOVA 
summary(test)[[1]][["Pr(>F)"]]

Post-hoc:

TukeyHSD(aov(y ~ x, data))

对于ANCOVA的post hoc:

TukeyHSD(aov(y ~ x+covariate, data),which='x')

如果ANOVA不满足正态性，用Kruskal-Wallis非参数检验：

kruskal.test(y~x, data = df) # post-hoc for Kruskal-Wallis Test： Steel-Dwass test library(FSA) dunnTest(y~x, data=df, method="bonferroni")

多重比较校正：

ps <- runif(20, max=.2) 
def <- p.adjust(ps) 
bon <- p.adjust(ps, method="bonferroni") 
bh <- p.adjust(ps, method="BH") #Benjamini, Hochberg (1995) default fdr 
by <- p.adjust(ps, method="BY") #Benjamini, Yekutieli (2001) another fdr 
x = matrix(rnorm(10000*5),nrow=10000) 
y = matrix(rnorm(10000*5),nrow=10000) 
p = sapply(1:10000, function(i) t.test(x[i,],y[i,])$p.val) 
q = p.adjust(p, method = "fdr")