R语言 数理统计

library(tidyverse)
library(readxl)

从excel文件的多个表里读入数据
excel_sheets()  #excel_sheets是readxl包的函数。该函数返回了porto_sheets.xlsx文件的所有sheet的名字

map_df(该文件各sheet的名称,read_excel, path = 文件名称)  #默认采用了按行纵向合并为一个数据框

合并函数map_df、dfr、dfc的时候,如果同一个文件的表合并,后面跟补充函数path


例子:
file <- "porto_sheets.xlsx"
file %>% excel_sheets() %>% 
  map_df(read_excel, path = file) -> porto_sheets 


从文件格式相同的多个文件里读入数据
library(fs)

dir_ls() 该函数返回了当前目录(工作目录)的文件清单

subdir <- "./porto"
filenames <- dir_ls(subdir, regexp = "\\.xlsx$")  读取该子目录的所有文件名
regexp: 其功能就是拼接文件路径及文件名字符串


一次性地就把该文件夹里的6个xlsx文件一次读入,并且纵向合并成一个数据框
filenames %>% map_df(read_excel) -> 
  	files

map_dfr(.x, .f): 返回数据框列表,再按行合并为一个数据框
map_dfc(.x, .f): 返回数据框列表,再 按列合并为一个数据框

library(pdftools)
mypdf <- "Wu Qian-Annals of Botany-Root-2018.pdf"
pdf_subset(mypdf, pages = 1, output = "aob1.pdf") 1页
pdf_subset(mypdf, pages = 5, output = "aob2.pdf") 5页
pdf_subset(mypdf, pages = 10:11, output = "aob3.pdf") 最后两页
执行后,在R的工作目录里能看到生成的3个pdf文件

可以计算pdf文件的页数
pdf_length("aob3.pdf")

可以将多个pdf文件重新合成一个新的pdf文件
pdf_combine(c("aob1.pdf","R_code.pdf", "aob3.pdf"), 
            output = "aob_new.pdf")

*******单因素方差分析——检验前提(正态性)

用W检验(H0: 样本数据与正态分布没有显著区别。H1: 样本数据与正态分布存在显著区别)
# 导入数据
df <- read.csv("sss.csv")   
# 进行正态性检验
shapiro.test(df[which(df[, 2] == "A"), ] $x)
shapiro.test(df[which(df[, 2] == "B"), ] $x)
shapiro.test(df[which(df[, 2] == "C"), ] $x)
shapiro.test(df[which(df[, 2] == "D"), ] $x)
shapiro.test(df[which(df[, 2] == "E"), ] $x)

W值越小,越接近 0,表示  样本数据越接近正态分布;
p值,如果 p-value小于显著性水平 α(0.05),则拒绝 H0


********单因素方差分析——检验前提(齐次性)

 (H0:各因子水平下的方差相同,H1:各因子水平下的方差不同)
bartlett.test(x~method, data = df)


********单因素方差分析——检验前提(离群点)
library(car)
outlierTest(aov(x~method, data=df))

方差分析模型可视为一种特殊的线性模型, 因此方差分析还可以使用第九章讲的线性模型函数lm( ), 并用函数anova( )提取其中的方差分析表, 因此aov(formula)等价于anova(lm(formula));

方差分析函数aov( )既适合于单因素方差分析, 也同样适用于双因素方差分析, 其中方差模型公式为x~A+B


**********多重比较
TukeyHSD(aov(x~method, data = df), ordered = F)
plot(TukeyHSD(aov(x~method, data = df)))


**********单因素协方差分析(ANCOVA)
协方差分析步骤如下:
(1)评估检验假设条件;(2)单因素协方差分析

aggregate函数可以按照要求把数据打组聚合,然后对聚合以后的数据进行加和、求平均等各种操作。
aggregate(litter$weight, by = list(litter$dose), FUN = mean)
by的功能类似于group_by
result:
>>aggregate(x[,3:4],by=list(sex=x$sex),FUN=mean)
  sex age height
1   F  26  152.5
2   M  30  168.0

library(HH)
ancova(weight~gesttime+dose, data = litter)

为了看协变量对组均值的影响,可以计算去除协变量效应后的组均值,可使用effects包中的effects()函数来计算调整的均值。

effect("dose", aov(weight~gesttime+dose, data = litter))

ANCOVA模型包含怀孕时间*剂量的交互项时,可对回归斜率的同质性进行检验,回归斜率同质性代码如下:
summary(aov(weight~gesttime*dose, data = litter))


***********双因素方差分析

R中的函数 aov():

不考虑交互作用x~A+B,加号表示两个因素具有可加性;
考虑交互作用的方差分析模型x~A+B+A:B(等价于A*B)
 
        ***无交互
# 导入数据
juice <- read.csv("juice.csv")
# 转化为因子型变量
juice$A <- factor(juice$A)
juice$B <- factor(juice$B)
# 双因素方差分析
juice.aov <- aov(X~A+B, data = juice)
summary(juice.aov)
检验方差齐次性
bartlett.test(X~A, data = juice) 
bartlett.test(X~B, data = juice) 
        ***有交互
# 导入数据   setwd(“D:/R-example”)
rats <- read.csv("rats.csv")
rats$Toxicant <- factor(rats$Toxicant)
rats$Cure <- factor(rats$Cure)

# 更改图形显示的界面排版
op <- par(mfrow = c(1, 2))
plot(Time~Toxicant+Cure, data = rats)
# 使用函数 interaction.plot()作出交互效应图, 考查因素之间交互作用是否存在:
with(rats, interaction.plot(Toxicant, Cure, Time, trace.label = "Cure"))
with(rats, interaction.plot(Cure, Toxicant, Time, trace.label = "Toxicant"))
曲线没有明显的相交情况出现,初步认为两个因素没有交互作用,方差分析函数 aov()进行因素间交互作用确认。
rats.aov <- aov(Time~Toxicant*Cure, data = rats)
summary(rats.aov)

进一步检验方差齐次性:
( bartlett和 Levene两种方法检验数据是否满足方差齐性的要求)
library(car)
leveneTest(rats$Time, rats$Toxicant)
leveneTest(rats$Time, rats$Cure)
bartlett.test(Time~Toxicant, data = rats)
bartlett.test(Time~Cure, data = rats)

**************线性回归

df_lb <- import('LA_biomass2.csv')
df_lb
#建立y关于x的回归方程
fitlm_lb <- lm(biomass~LA, data = df_lb)
summary(fitlm_lb)
#画散点图 abline加回归线
plot(biomass~LA, data = df_lb)
abline(fitlm_lb)



**************逐步回归

回归方程的系数的显著性不高, 有的甚至没有通过检验(X1与X2),这说明如果选择全部变量构造方程, 效果并不好. 这就涉及到变量选择的问题, 以建立“最优”的回归方程.
 变量选择与最优回归: R软件提供了获得“最优”回归方程的方法,“逐步回归法”的计算函数step( ), 它是以Akaike信息统计量为准则(简称AIC准 则), 通过选择最小的AIC信息统计量, 来达到删除或增加变量的目的. 

#导入数据
df_cd <- read.csv('ClimateData.csv')
#建立y关于x的回归方程
fitlm_cd <- lm(y~x1+x2+x3+x4+x5+x6, data = df_cd)  # +0
summary(fitlm_cd)
#用“一切子集回归法”来进行逐步回归
steplm_cd <- step(fitlm_cd, direction = 'both') 
summary(steplm_cd)
*Direction确定逐步搜索的方向,默认为both“一切子集回归法”,backward
是“向后法”,forward是“向前法”。

******************多元线性回归的假设前提
确定了回归模型的自变量并初步得到一个线性回归模型,并不是直接可以拿来用的,还要进行验证和诊断。

#残差分析和异常点检测:
residuals(): 普通残差
rstandard(): 标准化残差
rstudent():学生化残差

# 一般把标准化残差的绝对值 >=2的观测值看作可疑点,把标准化残差的绝对值 >=3的观测值看作异常点

#计算残差
y_res <- residuals(steplm_cd)
#计算标准化
y_rst <- rstandard(steplm_cd)
print(y_rst)
#计算预测值
y_fit <- predict(steplm_cd)
#画残差散点图
plot(y_res~y_fit)
plot(y_rst~y_fit)

#方差稳定变换(对数变换)
steplm_cd_up <- update(steplm_cd, log(.)~.)

y_rst_up <- rstandard(steplm_cd_ up)
y_fit_up <- predict(steplm_cd_up)
plot(y_rst_up ~y_fit_up)

#剔除一些点再进行残差分析
fitlm_cd_new <- lm(log(y)~x1+x2+x3+x4+x5+x6, data = df_cd[-c(3,12),] )
steplm_cd_new <- step(fitlm_cd_new)
y_rst_new <- rstandard(steplm_cd_new)
y_fit_new <- predict(  steplm_cd_new)
plot(y_rst_new~y_fit_new)

#函数 plot()和 influence.measures()可以用来绘制诊断图,计算诊断统计量。

par(mfrow = c(2,2))
plot(steplm_cd)
influence.measures(steplm_cd)

*******回归预测
#回归预测分为点预测和区间预测,R中用 predict()实现
preds <- data.frame(x1 = 29, x2 = 1,x3 = 9, x4 = 2, x5 = 6, x6 = 265)
predict(steplm_cd, newdata = preds, interval = 'prediction', level = 0.95)


car包里的influencePlot()函数能一次性同时检查离群点、高杠杆点、强影响点。
#导入数据并建立y关于x的逐步回归方程
library(car)
df_cd <- read.csv("ClimateData.csv")
fitlm_cd <- lm(y~x1+x2+x3+x4+x5+x6, data = df_cd)
steplm_cd <- step(fitlm_cd) 
#绘制回归影响图
influencePlot(fitlm_cd,id.method = "identity", main = "Influence Plot", sub = "Circle size is proportional to Cook's distance")
多重共线性 
    理想中的线性模型各个自变量应该是线性无关的,若自变量间有些相关,就有可能出现多重共线性。
    如果共线性严重时,模型或数据的微小变化有可能造成系数估计的较大变化,则模型不稳定。
    一般用方差膨胀因子 VIF (Variance Inflation Factor)来衡量共线性。一般认为 VIF超过 5或 10,则有多重共线性问题。另一种是条件数 (condition number),用 k表示。一般 k >15时有共线性问题,k >30时说明共线性问题严重。理想中的线性模型 VIF = 1,表示完全不存在共线性。

rm(list = ls())
library(tidyverse)
library(rio)
library(nycflights13)
library(ggplot2)
library(lubridate)
library(readxl)
library(fs)
library(multcomp)
library(HH)
library(car)

df_cd <- read.csv("ClimateData.csv")
fitlm_cd <- lm(y~x1+x2+x3+x4+x5+x6, data = df_cd)
steplm_cd <- step(fitlm_cd) 
#检验多重共线性
vif(fitlm_cd )
vif(steplm_cd)


#######查看自变量间的相关系数
plot(longley[, 2:7])
cor(longley[, 2:7])

#######用岭回归处理多重共线性问题
library(ridge)
a  <- linearRidge(GNP.deflator ~ .,data = longley)
summary(a)

#######用偏最小二乘回归处理多重共线性问题
library(pls)
ap <- plsr(GNP.deflator ~ ., data = longley)
summary(ap)

#root mean squared error of prediction
RMSEP(ap)
#the mean squared error of prediction
MSEP(ap)
R2(ap)
coef(ap)

#将RMSE、MSE、R²图并列展示par(mfrow = c(1,3))
plot(RMSEP(ap)) 
plot(MSEP(ap)) 
plot(R2(ap))

#主成分分析

(1)princomp(x, cor=F, scores=T)
其中, x为主成分分析数据集,cor=T和F分别代表是基于相关系数矩阵计算还是样本协方差矩阵(数据标准化)计算。scores则代表是否存储主成分得分。

(2) FactoMineR包的PCA() 
PCA(X, scale.unit = T, ncp =, graph = T)
其中,x为主成分分析数据集, scale.unit=T代表是分析前标准化数据。ncp为主成分个数, graph = T代表显示图像。

(3) factoextra包
get_eigenvalue(): 提取主成分的特征值/方差
get_pca_ind(): 提取个体结果
get_pca_var(): 提取变量结果


可视化主成分分析结果函数:
1)screeplot():可视化主成分特征值,碎石图
2)factoextra包
fviz_eig(): 可视化主成分特征值,碎石图——常用
fviz_pca_ind():可视化个体结果——常用
fviz_pca_var():可视化变量结果
fviz_pca_biplot(): 主成分分析分类作图

第一种方法:
#删除因子型变量
PCA <- princomp(iris[,-5])
summary(PCA)
#碎石图
screeplot(PCA, type="line")

第二种方法:
library(FactoMineR)
iris_pca <- PCA(iris[,-5], graph =F)
#提取特征值
library(factoextra)
eig_val <- get_eigenvalue(iris_pca)
eig_val
#碎石图
scree_plot <-fviz_eig(iris_pca)
scree_plot


#个体PCA可视化
ind_plot <- fviz_pca_ind(iris_pca,
        geom.ind = "point", # 点图
        col.ind = iris$Species, # 颜色分类
        palette = c("#00AFBB", "#E7B800", "#FC4E07"),#颜色设置
        addEllipses = TRUE, # 添加椭圆
        legend.title = "Groups") #添加标题
ind_plot

聚类分析11个典型的步骤:
1. 选择合适的变量。选择你感觉可能对识别和理解数据中不同观测值分组有重要影响的变量。
2. 缩放数据。标准化数据,最常用的方法是将每个变量标准化为均值0和标准差为1的变量,常用scale()函数。
3. 寻找异常点。许多聚类方法对异常值十分敏感。
·可以通过outliers包中的函数来筛选异常单变量离群点。
·mvoutlier包能识别多元变量的离群点的函数。
·使用对异常值稳健的聚类方法划分聚类。

4.计算距离
     虽然所使用的算法差异大,但是通常都需要计算被聚类的实体之间的距离。最常用欧几里得距离,其他可选曼哈顿距离、兰式距离、非对称二元距离、最大距离和闵可夫斯基距离。
      R中自带的dist()函数能够用来计算矩阵或数据框中所有行之间的距离。格式是:
                                      dist(x, method=)
x表示输入数据,并且默认为欧几里得距离。函数默认返回一个下三角矩阵,但是as.matirx()函数可使用标准括号符号得到距离。

5.选择聚类算法(层次聚类、快速聚类、有序聚类)。
6. 获得一种或多种聚类方法。使用步骤(5)选择的方法。
7. 确定类的数目。NbClust包提供了众多的指数来确定一个聚类分析里类的最佳数目。不能保证这些指标得出的结果都一致,但是可以作为选择聚类个数K值的一个。
8. 获得最终的聚类解决方案
9. 结果可视化
10. 解读类。
11. 验证结果。采用不同的聚类方法是否会产生相同的类?fcp, clv, clValid包包含了评估聚类解的稳定性的函数

#层次聚类
library('flexclust')
data(nutrient)
row.names(nutrient) <- tolower(row.names(nutrient))
head(nutrient)
#由于能量变化范围比其他变量更大,缩放数据有利于均衡各变量的影响。
nutrient_scaled <- scale(nutrient)
#用来计算矩阵或数据框中所有行之间的距离,默认欧几里得距离
d <- dist(nutrient_scaled) #默认下三角矩阵
#先不考虑分类数,层次聚类
fit_average <- hclust(d, method="average")
rect.hclust (fit, k=) #设定分类数用这个函数
#hang是表明谱系图中各类所在的位置,hang取负值时,类从底部画起
plot(fit_average, hang = -1,  main = "Average Linkage Clustering")

#快速聚类
当样本容量n较大时,用层次聚类法聚类有几点不足:
(1)距离矩阵的计算量很大;
(2)样品一旦归到某个类就不能改变;
(3)不同的聚类方法,产生的聚类结果可能不同。


算法:
1.选择K个中心点(随机选择K行);
2.把每个数据点分配到离它最近的中心点;
3.重新计算每类中的点到该类中心点距离的平均值;
4.分配每个数据到它最近的中心点;
5.重复步骤3,4直到所有观测值不再被分配或是达到最大的迭代次数(R默认10次)。

代码:
#确定分类数
k_means <- kmeans(d,3)  # 可用pam(data, k, metric="euclidean")代替
#其中,data为数据,k为分类数, metric表示使用的相似性/相异性的度量,默认欧几里得距离 。

#可视化
install.packages(' factoextra ')
library(factoextra)
fviz_cluster(k_means, nutrient)
#判别分析

判别方法:1)距离判别法 2)Fisher线性判别法 3)Bayes判别法

#Fisher线性判别法:
library('MASS')
data(iris)
head(iris)
colnames(iris) <- tolower(colnames(iris))
head(iris)
#Fisher判别分析
ld <- lda(species~sepal.length+sepal.width+petal.length+petal.width,data=iris)
#对原始数据进行回判分类
p_iris <- predict(ld) 
new_class<- p_iris$class
# p_iris$x 给出了p_iris中两个判别函数对应的值
cbind(iris$species, new_class, p_iris$x) 

#Bayes判别法

  • 2
    点赞
  • 20
    收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:1024 设计师:我叫白小胖 返回首页
评论

打赏作者

我叫杨傲天

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值