本文基于R语言进行基本数据统计分析,包括基本作图,线性拟合,逻辑回归,bootstrap采样和Anova方差分析的实现及应用。
不多说,直接上代码,代码中有注释。
1. 基本作图(盒图,qq图)
#basic plot
boxplot(x)
qqplot(x,y)
2. 线性拟合
#linear regression
n = 10
x1 = rnorm(n)#variable 1
x2 = rnorm(n)#variable 2
y = rnorm(n)*3
mod = lm(y~x1+x2)
model.matrix(mod) #erect the matrix of mod
plot(mod) #plot residual and fitted of the solution, Q-Q plot and cook distance
summary(mod) #get the statistic information of the model
hatvalues(mod) #very important, for abnormal sample detection
3. 逻辑回归
#logistic regression
x <- c(0, 1, 2, 3, 4, 5)
y <- c(0, 9, 21, 47, 60, 63) # the number of successes
n <- 70 #the number of trails
z <- n - y #the number of failures
b <- cbind(y, z) # column bind
fitx <- glm(b~x,family = binomial) # a particular type of generalized linear model
print(fitx)
plot(x,y,xlim=c(0,5),ylim=c(0,65)) #plot the points (x,y)
beta0 <- fitx$coef[1]
beta1 <- fitx$coef[2]
fn <- function(x) n*exp(beta0+beta1*x)/(1+exp(beta0+beta1*x))
par(new=T)
curve(fn,0,5,ylim=c(0,60)) # plot the logistic regression curve
3. Bootstrap采样
# bootstrap
# Application: 随机采样,获取最大eigenvalue占所有eigenvalue和之比,并画图显示distribution
dat = matrix(rnorm(100*5),100,5)
no.samples = 200 #sample 200 times
# theta = matrix(rep(0,no.samples*5),no.samples,5)
theta =rep(0,no.samples*5)
for (i in 1:no.samples)
{
j = sample(1:100,100,replace = TRUE)#get 100 samples each time
datrnd = dat[j,]#select one row each time
lambda = princomp(datrnd)$sdev^2#get eigenvalues
# theta[i,] = lambda
theta[i] = lambda[1]/sum(lambda)#plot the ratio of the biggest eigenvalue
}
# hist(theta[1,]) #plot the histogram of the first(biggest) eigenvalue
hist(theta)#plot the percentage distribution of the biggest eigenvalue
sd(theta)#standard deviation of theta
#上面注释掉的语句,可以全部去掉注释并将其下一条语句注释掉,完成画最大eigenvalue分布的功能
4. ANOVA方差分析
#Application:判断一个自变量是否有影响 (假设我们喂3种维他命给3头猪,想看喂维他命有没有用)
#
y = rnorm(9)#weight gain by pig(Yij, i is the treatment, j is the pig_id), 一般由用户自行输入
#y = matrix(c(1,10,1,2,10,2,1,9,1),9,1)
Treatment <- factor(c(1,2,3,1,2,3,1,2,3)) #each {1,2,3} is a group
mod = lm(y~Treatment) #linear regression
print(anova(mod))
#解释:Df(degree of freedom)
#Sum Sq: deviance (within groups, and residuals) 总偏差和
# Mean Sq: variance (within groups, and residuals) 平均方差和
# compare the contribution given by Treatment and Residual
#F value: Mean Sq(Treatment)/Mean Sq(Residuals)
#Pr(>F): p-value. 根据p-value决定是否接受Hypothesis H0:多个样本总体均数相等(检验水准为0.05)
qqnorm(mod$residual) #plot the residual approximated by mod
#如果qqnorm of residual像一条直线,说明residual符合正态分布,也就是说Treatment带来的contribution很小,也就是说Treatment无法带来收益(多喂维他命少喂维他命没区别)
如下面两图分别是
(左)用 y = matrix(c(1,10,1,2,10,2,1,9,1),9,1)和
(右)y = rnorm(9)
的结果。可见如果给定猪吃维他命2后体重特别突出的数据结果后,qq图种residual不在是一条直线,换句话说residual不再符合正态分布,i.e., 维他命对猪的体重有影响。
@[toc]
假设检验的前提是要满足正态分布和方差齐性
组内平方和SSE:同一组内的数据误差平方和
组间平方和SSA:不同组之间的数据误差平方和
一个分类型自变量
例如四个班级学生的语文成绩,班级是分类型自变量,四个班级是自变量的四个水平
测试班级对成绩的影响
因为p<0.001,说明班级对成绩的影响非常显著
图中跨越0分界线的班级对,有较大可能落在0上,也就是说两个班级之间没有明显差异。其他班级说明都有明显差异。
同一班级在大学三年的三次测试
p<0.001,说明学生成绩在大学三年中有显著差异。球形检验的p-value大于0.05,所以可以认为方差相等。
Mauchly's Test for Sphericity :适用于重复测量时检验不同测量之间的差值的方差是否相等,用于三次以及三次之上。
Sphericity Corrections :球形矫正,当方差不相等时进行矫正,矫正方法有the Greenhouse-Geisser (1959), the Huynh-Feldt (1976), 简称GG和HF。
两个分类型自变量
例如探究 词汇量 和 话题熟悉度 对学生作文成绩的影响
词汇量和话题熟悉度两个变量对成绩的影响都很显著,交互项对成绩影响不显著。
探究班级和测试次数对学生成绩的影响
班级和测试次数在原始检验中都很显著,然后交叉项不显著。
但是在球形检验中,推翻了方差齐性的假设,所以tests需要使用球形矫正之后的p值,classes不用。
矫正之前tests的p-value = 3.482406e-04,矫正之后的p-value = 0.001左右。
R语言与统计-1:t检验与秩和检验
方差分析适用于多组均数的比较 (在完全随机设计的实验中,两组均数的t检验和方差分析是完全等价的。但t检验只能用于两组的均数比较,对于三组和三组以上的均数比较,就需要用到方差分析。)
可以看到这个数据集只有两个变量,其中治疗是分类变量(因子型),有5个水平。response是数值型变量。要对每种治疗所对应的response的均值进行比较,就只能用方差分析而不能用t检验。
符合正态分布
要比较均值的数据写~左边,分组变量写右边。p=0.9653,方差齐。
写法同上,方差齐。
需要注意的是,如果检验出方差不齐,我们第一步不是立马选择进行非参数检验,而是首先要判断有无异常值存在,因为异常值对方差的影响很大。当然,到这一步才来检验有无异常值是不符合数据分析的流程的,异常值在进行数据初步处理的时候就因该被发现和处理掉。
方差分析包括 单因素方差分析 , 多因素方差分析 , 协方差分析 , 多元方差分析 , 重复测量数据方差分析 。
gplots包的plotmeans函数 对上述结果进行可视化
使用ToothGrowth数据集进行演示
aov函数
不考虑supp和dose之间的交互作用的情况。结果显示两个因素都对小鼠牙齿生长影响显著。
考虑两个因素之间的交互作用:将上面的+换成*。结果显示两个因素都对小鼠牙齿生长影响显著而且两者间的互相影响也不容忽视。
可视化
上述结果已经知道了再五组数据中的均值不全相等,下一步想知道哪些相等哪些不等,就要对这五组进行两两比较。
输出的结果:从左往右依次是:两两比较、两两间的差值、lwr是95%可信线的下限,upr是上限。最后是p值。
将结果可视化:
线段中点是均值,两端是95%置信区间,跨过0说明没有显著差异。
在进行方差分析时,所有混杂因素统称为 协变量 。
检验dose对weight的影响。出生时间gesttime是协变量。
aov后面小括号里写的顺序:结果变量~协变量+自变量。如果要看协变量和自变量之间是否存在交互,在后面写+协变量:自变量。最后是data=数据集。
结果显示两个变量之间不存在交互效应(p=0.17889, >0.05),可以认为它们的斜率是相同的。
因变量不止一个,但是需要将它们作为一个整体同时进行分析。例如:某种药物对患者血红蛋白浓度,红细胞计数,外周血细胞因子水平等多种因素的影响。
使用 manovs()函数 进行性多元方差分析
参考: https://blog.csdn.net/dingming001/article/details/72822270