【R语言入门与数据分析-5】 数据分析实战

Python017

【R语言入门与数据分析-5】 数据分析实战,第1张

老师的吐槽大会,乐死我了。hhh

regression,通常指用一个或者多个预测变量,也称自变量或者解释变量,来预测响应变量,也称为因变量、效标变量或者结果变量的方法

存在多个变量

AIC 考虑模型统计拟合度、用来拟合的参数数目

AIC值越小,越好

更多的变量:

图一:是否呈线性关系, 是

图二:是否呈正态分布,一条直线,正态分布

图三:位置与尺寸图,描述同方差性,如果方差不变,水平线周围的点应该是随机分布

图四:残差与杠杆图,对单个数据值的观测,鉴别离群点、高杠杆点、强影响点

模型建好,用predict函数对剩余500个样本进行预测,比较残差值,若预测准确,说明模型可以。

analysis of variance,简称ANOVA,也称为变异数分析。用于两个及两个以上样本均数差别的显著性检验。广义上,方差分析也是回归分析的一种,只不过线性回归的因变量一般是连续型变量。自变量是因子时,研究关注的重点通常会从预测转向不同组之间的差异比较。也就是方差分析。

power analysis,可以帮助在给定置信度的情况下,判断检测到给定效应值所需的样本量。也可以在给定置信度水平情况下,计算在某样本量内能检测到给定效应值的概率

拓展了线性模型的框架,包含了非正态因变量的分析。线性回归、方差分析都是基于正态分布的假设

-泊松回归 ,用来为计数资料和列联表建模的一种回归分析。泊松回归假设因变量是泊松分布,并假设它平均值的对数可被未知参数的线性组合建模

-logistic 回归

通过一系列连续型或者类别型预测变量来预测二值型结果变量是,logistic 回归是一个非常有用的工具。流行病学研究中用的多。

Principal Component Analysis,PCA,探索和简化多变量复杂关系的常用方法。 是一种数据降维技巧。可以将大量相关变量转化为一组很少的不相关变量。这些无关变量成为主成分。主成分是对原始变量重新进行线性组合,将原先众多具有一定相关性的指标,重新组合为一组的心得相互独立的综合指标。

探索性因子分析法 exploratory factor analysis,简称为EFA,是一系列用来发现一组变量的潜在结构的方法。通过找寻一组更小的、潜在的活隐藏的结构来解释已观测到的、显式的变量间的关系

因子分析步骤与PCA一致

啤酒与尿布

R语言中的多元方差分析

1、当因变量(结果变量)不止一个时,可用多元方差分析(MANOVA)对它们同时进行分析。

library(MASS)

attach(UScereal)

y <- cbind(calories, fat, sugars)

aggregate(y, by = list(shelf), FUN = mean)

Group.1 calories fatsugars

1 1 119.4774 0.6621338 6.295493

2 2 129.8162 1.3413488 12.507670

3 3 180.1466 1.9449071 10.856821

cov(y)

calories fat sugars

calories 3895.24210 60.674383 180.380317

fat60.67438 2.713399 3.995474

sugars180.38032 3.995474 34.050018

fit <- manova(y ~ shelf)

summary(fit)

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

shelf 1 0.195944.955 3 61 0.00383 **

Residuals 63

---

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

summary.aov(fit)

Response calories :

Df Sum Sq Mean Sq F valuePr(>F)

shelf1 45313 45313 13.995 0.0003983 ***

Residuals 63 2039823238

---

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

Response fat :

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

shelf1 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)

shelf1 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)多元正态性

第一个假设即指因变量组合成的向量服从一个多元正态分布。可以用Q-Q图来检验该假设条件。

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 = "QQ

Plot Assessing Multivariate Normality",

ylab = "Mahalanobis D2")

abline(a = 0, b = 1)

identify(coord$x, coord$y, labels = row.names(UScereal))

如果所有的点都在直线上,则满足多元正太性。

2、方差—协方差矩阵同质性即指各组的协方差矩阵相同,通常可用Box’s M检验来评估该假设

3、检测多元离群点

library(mvoutlier)

outliers <- aq.plot(y)

outliers

数据准备

许多实际情况中统计假设(假定观测数据抽样自正态分布或者其他性质较好的理论分布)并不一定满足,比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况,这时基于随机化和重抽样的统计方法就可派上用场。

置换检验的定义

置换检验(Permutation test),也称随机化检验或重随机化检验,是Fisher于20世纪30年代提出的一种基于大量计算(computationally intensive),利用样本数据的全(或随机)排列,进行统计推断的方法,因其对总体分布自由,应用较为广泛,特别适用于总体分布未知的小样本资料,以及某些难以用常规方法分析资料的假设检验问题。

置换检验的原理

1、提出原假设,比如XX处理后结果没有变化

2、计算统计量,如两组的均值之差,记作t0

3、将所有样本放在一起,然后随机排序进行分组,再计算其统计量t1

4、重复第3步骤,直至所有排序可能性都齐全(比如有A组有n样本,B组有m样本,则总重复次数相当于从n+m中随机抽取n个的次数),得到一系列的统计量(t1-tn)

5、最后将这些统计量按照从小到大排序,构成抽样分布,再看t0是否落在分布的置信区间内(如95%置信区间),这时候可计算一个P值(如果抽样总体1000次统计量中大于t0的有10个,则估计的P值为10/1000=0.01),落在置信区间外则拒绝原假设

6、如果第3步骤是将所有可能性都计算了的话,则是精确检验;如果只取了计算了部分组合,则是近似结果,这时一般用蒙特卡罗模拟(Monte Carlo simulation)的方法进行置换检验

7、置换检验和参数检验都计算了统计量,但是前者是跟置换观测数据后获得的经验分布进行比较,后者则是跟理论分布进行比较。

请牢记:置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样(当做近似检验时)。因此,每次检验的结果都有所不同。

coin包提供了一个进行置换检验的一般性框架。通过该包,你可以回答如下问题。

响应值与组的分配独立吗?

两个数值变量独立吗?

两个类别型变量独立吗?

表12-2列出来的每个函数都是如下形式:

function_name(formula, data, distribution=)

其中:

 formula描述的是要检验变量间的关系。示例可参见表12-2;

 data是一个数据框;

 distribution指定经验分布在零假设条件下的形式,可能值有exact,asymptotic和

approximate。

若distribution = "exact",那么在零假设条件下,分布的计算是精确的(即依据所有可能的排列组合)。当然,也可以根据它的渐进分布(distribution = "asymptotic")或蒙特卡洛重抽样(distribution = "approxiamate(B = #)")来做近似计算,其中#指所需重复的次数。

distribution = "exact"当前仅可用于两样本问题。

传统t检验表明存在显著性差异(p <0.05),而精确检验却表明差异并不显著(p >0.072)。

第7章我用自己的数据进行了t检验,对比一下传统t检验和置换检验,结果如下:

两种检验方式下结果都是显著的

Wilcoxon-Mann-Whitney U检验

coin包规定所有的类别型变量都必须以因子形式编码。

wilcox.test()默认计算的也是精确分布。

K样本检验的置换检验

通过chisq_test()或cmh_test()函数,我们可用置换检验判断两类别型变量的独立性。 当数据可根据第三个类别型变量进行分层时,需要使用后一个函数。若变量都是有序型,可使用lbl_test()函数来检验是否存在线性趋势。

卡方独立性检验

卡方独立性检验的置换检验

你可能会有疑问,为什么需要把变量Improved从一个有序因子变成一个分类因子?(好问题!)这是因为,如果你用有序因子,coin()将会生成一个线性与线性趋势检验,而不是卡方检验。

结果解读:两种检验下p值都是小于0.05,说明Treatment和Improved之间相互不独立

自己数据的演示

结果解读:p值均为1,表明nitrogen和variety相互独立。

spearman_test()函数提供了两数值变量的独立性置换检验。

当处于不同组的观测已经被分配得当,或者使用了重复测量时,样本相关检验便可派上用场。

对于两配对组的置换检验,可使用wilcoxsign_test()函数;多于两组时,使用friedman_test()函数。

自己数据演示

lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版,能够进行置换检验,而非正态理论检验。

lmp()和aovp()函数的参数与lm()和aov()函数类似,只额外添加了perm =参数。

perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样,直至估计的标准差在估计的p值0.1之下,判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意,若观测数大于10,perm = "Exact"将自动默认转为perm = "Prob",因为精确检验只适用于小样本问题。

简单线性回归的置换检验

R语言实战的例子:

多项式回归的置换检验

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

当两种方法所得结果不一致时,你需要更加谨慎地审视数据,这很可能是因为违反了正态性假设或者存在离群点。

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

值得注意的是,当将aovp()应用到方差分析设计中时,它默认使用唯一平方和法(SAS也称为类型III平方和)。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和(SAS是类型I平方和)。每种效应依据模型中先出现的效应做相应调整。对于平衡设计,两种方法结果相同,但是对于每个单元格观测数不同的不平衡设计,两种方法结果则不同。不平衡性越大,结果分歧越大。若在aovp()函数中设定seqs = TRUE,可以生成你想要的序贯平方和。

你可能已经注意到,基于正态理论的检验与上面置换检验的结果非常接近。在这些问题中数据表现非常好,两种方法结果的一致性也验证了正态理论方法适用于上述示例。

当然,置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大)、存在离群点、样本很小或无法做参数检验等情况。不过,如果初始样本对感兴趣的总体情况代表性很差,即使是置换检验也无法提高推断效果。

置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。不过,置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是,这正是自助法大显神通的地方。

所谓自助法,即从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。 无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设。

倘若你假设均值的样本分布不是正态分布,该怎么办呢?可使用自助法。

(1)从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一直都不会被选中。

(2)计算并记录样本均值。

(3)重复1和2一千次。

(4)将1000个样本均值从小到大排序。

(5)找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限定了95%的置信区间。

boot包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。

一般来说,自助法有三个主要步骤。

(1)写一个能返回待研究统计量值的函数。如果只有单个统计量(如中位数),函数应该返回一个数值;如果有一列统计量(如一列回归系数),函数应该返回一个向量。

(2)为生成R中自助法所需的有效统计量重复数,使用boot()函数对上面所写的函数进行处理。

(3)使用boot.ci()函数获取第(2)步生成的统计量的置信区间。

主要的自助法函数是boot(),它的格式为:

bootobject <- boot(data=, statistic=, R=, ...)

参数见下表:

boot()函数调用统计量函数R次,每次都从整数1:nrow(data)中生成一列有放回的随机指标,这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算,结果存储在bootobject中。

你可以用bootobject t0和bootobject t来获取这些元素。

一旦生成了自助样本,可通过print()和plot()来检查结果。如果结果看起来还算合理, 使用boot.ci()函数获取统计量的置信区间。格式如下:

boot.ci(bootobject, conf=, type= )

type参数设定了获取置信区间的方法。perc方法(分位数)展示的是样本均值,bca将根据偏差对区间做简单调整。

回归的R平方值

1000次自助抽样

输出结果

结果可视化

95%的置信区间获取

回归系数向量函数

自助抽样1000次

获得车重和发动机排量95%的置信区间

置换检验和自助法并不是万能的,它们无法将烂数据转化为好数据。当初始样本对于总体情况的代表性不佳,或者样本量过小而无法准确地反映总体情况,这些方法也是爱莫能助。

参考资料: