pvrect r语言 聚类_R语言回归分析及相关分析
#两样本问题t检验、方差齐次性检验、分布检验
#相关与回归分析
#例一某种元件的寿命X(小时)，服从正态分布，N(mu,sigma2)，其中mu,sigma2均未知，16只元件的寿命如下：问是否有理由认为元件的平均寿命大于225小时。
#分析：按题意，需检验
#H0： μ ≤ 225 H1: μ > 225
#此问题属于单边检验问题，可以使用R语言t.test
#t.test(x,y=NULL,
alternative=c(“two.sided”,“less”,“greater”),
mu=0,paired=FALSE,var.equal=FALSE,
conf.level=0.95)
#其中x,y构成两向量，(如果只提供x,则作单个正态总体的均值检验，如果提供x,y则作两个总体的均值检验)，
#alternative表示备择假设，two.sided(缺省)，双边检验(H1:μ≠H0),less表示单边检验(H1:μμ0)，
#mu表示原假设μ0，paired表示是否配对样本T检验，conf.level置信水平，即1-α，通常是0.95，var.equal是逻辑变量，var.equal=TRUE表示两样品方差相同，var.equal=FALSE(缺省)表示两样本方差不同。
X
t.test(X, alternative = “greater”, mu = 225)
#可见P值为0.257 > 0.05 ，不能拒绝原假设，接受H0，即平均寿命不大于225小时。
#例二 在平炉上进行的一项试验以确定改变操作方法的建议是否会增加钢的得率，试验时在同一个平炉上进行的，每炼一炉刚时除操作方法外，其它条件都尽可能做到相同，先用标准方法炼一炉，然后用新方法炼一炉，以后交替进行，各炼了10炉，其得率分别为
#标准方法 78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3
#新方法 79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1
#设这两个样本相互独立，且分别来自正态总体N(μ1, σ2)和N(μ2, σ2)，其中μ1，μ2和σ2未知。问新的操作能否提高得率？(取α=0.05)
#分析：根据题意，需要假设
#H0： μ1 ≥ μ2 H1: μ1 < μ2
#这里假定σ12=σ22=σ2，因此选择t.test，var.equal=TRUE
#R代码：
X
Y
t.test(X,Y,var.equal = TRUE,alternative = “less”)
#可见P值<0.05，拒绝原假设，即新的操作能够提高得率。
#因为数据是成对出现的，所以采用成对数据t检验比上述的双样本均值检验更准确。
#所谓成对t检验就是Zi=Xi-Yi，再对Z进行单样本均值检验(与0比较)
X
Y
t.test(X-Y, alternative = “less”)
t.test(X,Y,paired = TRUE,alternative = “less”)
#例二方差相同吗？检验方差是否相同
#H0： σ12 = σ22 H1: σ12 ≠ σ22
#方差检验可以用var.test
#var.test(x, y, ratio = 1,
alternative = c(“two.sided”, “less”, “greater”),
conf.level = 0.95, …)
#x,y是来自两样本数据构成的向量，ratio是方差比的原假设，缺省值为1.alternative是备择假设，two.sided表示双边检验,greater表示单边检验
#R代码：
X
Y
var.test(X,Y)
#可见P值为0.559>0.05，不拒绝原假设，认为两者方差相同
#二项分布的参数检验
#有一批蔬菜种子的平均发芽率p0=0.85，现随即抽取500粒，用种衣剂进行浸种处理，结果有445粒发芽。试检验种衣剂对种子发芽率有无效果。
#根据题意，所检验的问题为
#H0：p=p0=0.85， H1：p≠p0
#可以用R语言的binom.test
#binom.test(x, n, p = 0.5,
alternative = c(“two.sided”, “less”, “greater”),
conf.level = 0.95)
#其中x是成功的次数；或是一个由成功数和失败数组成的二维向量。n是试验总数，当x是二维向量时，此值无效。P是原假设的概率。
#R语言代码：
binom.test(445,500,p=0.85)
#可知P值0.01207<0.05，拒绝原假设，说明种衣剂对种子的发芽率有显著效果。
#拟合优度检验
#例：某体育网站为了确定球迷对5支具备夺取世界杯冠军球队的支持情况，随即抽取了1000名球迷进行问卷调查：每个人选一支最支持的队伍。这5支球队分别记为A、B、C、D、E。下表是根据样本资料整理的5支球队支持人数的频数分布。
#试根据这些数据判断球迷对五支球队的支持数量有无明显差异？
最喜欢的球队 A B C D E
支持的人数 210 312 170 85 223
#分析：
#如果球迷对5支球队的支持无显著差异，那么，就可以认为喜好这5支球队的人呈均匀分布，即支持每队的人数各占20%。据此假设
#H0：支持5支球队的人数服从分布均匀
#可以使用Pearson χ2拟合优度检验，R语言中调用chisq.test(X)
#chisq.test(x, y = NULL, correct = TRUE,
p = rep(1/length(x), length(x)), rescale.p = FALSE,
simulate.p.value = FALSE, B = 2000)
#其中x是由观测数据构成的向量或者矩阵，y是数据向量(当x为矩阵时，y无效)。correct是逻辑变量，标明是否用于连续修正，TRUE(缺省值)表示修正，FALSE表示不修正。p是原假设落在小区间的理论概率，缺省值表示均匀分布，rescale.p是逻辑变量，选择FALSE(缺省值)时，要求输入的p满足和等于1；选择TRUE时，并不要求这一点，程序将重新计算p值。simulate.p.value逻辑变量(缺省值为FALSE)，当为TRUE，将用仿真的方法计算p值，此时,B表示仿真的此值。
#R语言代码：
X
chisq.test(X)
#p值非常接近于0，拒绝原假设
#泊松分布检验
#为研究电话总机在某段时间内接到的呼叫次数是否服从Poisson分布，现收集了42个数据，如下表所示，通过对数据的分析，
#问能否确认在某段时间内接到的呼叫次数服从Poisson分布(α = 0.1)？
#接到呼叫次数 0 1 2 3 4 5 6
#频数 7 10 12 8 3 2 0
X
#计算理论分布，其中mean(rep(X,Y))为样本均值
q
p=rep(0,n)
p[1]
for (i in 2:(n-1))
p[i]
#作检验
chisq.test(Y,p=p)
#提示结果可能不准确，因为皮尔森卡方拟合由度检验要求分组后每组的频数至少要大于等于5，而后三组中出现的频率分别为3，2，0，均小于5，
#解决问题的方法是将后三组合成一组，此时的频数为5，满足要求，重写R语言代码
#输入数据
X
#计算理论分布，其中mean(rep(X,Y))为样本均值
q
p
p[1]
for (i in 2:(n-1))
p[i]
#重新分组
Z
#重新计算理论分布
n
#作检验
chisq.test(Z,p=p)
#p值大于0.1，不拒绝原假设，表明服从泊松分布
#正态分布的检验
#已知15名学生体重如下，问是否服从正态分布
#75.0, 64.0, 47.4, 66.9, 62.2, 62.2, 58.7, 63.5,
#66.6, 64.0, 57.0, 69.0, 56.9, 50.0, 72.0
w
66.6, 64.0, 57.0, 69.0, 56.9, 50.0, 72.0)
shapiro.test(w)
#P值>0.05，接受原假设，认为来自正态分布总体。
#也可以用ks.test检验正态性
#统计学里, Kolmogorov–Smirnov 检验(亦称：K–S 检验)是用来检验数据是否符合某种分布的一种非参数检验，通过比较一个频率分布f(x)与理论分布g(x)或者两个观测值分布来判断是否符合检验假设。其原假设H0:两个数据分布一致或者数据符合理论分布。拒绝域构造为：D=max| f(x)- g(x)|，当实际观测值D>D(n,α)则拒绝H0，否则则接受H0假设。由于KS检验不需要知道数据的分布情况，在小样本的统计分析中效果比较好。(大样本数据下，使用t-检验；小样本数据，使用t-检验会出现较大的偏差)
#ks.test(x, y, …,
alternative = c(“two.sided”, “less”, “greater”),
exact = NULL)
#x:观测值向量
#y:第二观测值向量或者累计分布函数，如pnorm(正态分布函数，一般做正态检测的时候直接输入pnorm)，只对连续CDF有效
#alternative = c(“two.sided”, “less”, “greater”):双侧检验还是单侧检验
#exact：默认为NULL，也可以是其他逻辑值，表明是否需要计算精确的P值
#结果解释：D：D值越小，越接近0，表示样本数据越接近正态分布(简单来说，D越小越好),p：p-value小于显著性水平α(0.05)，则拒绝H0(p越大越好)
x1=c(87,77,92,68,80,78,84,77,81,80,80,77,92,86,76,80,81,75,77,72,81,90,84,86,80,68,77,87,76,77,78,92,75,80,78)
x1
ks.test(x1,pnorm)
#注意事项：
#在做单样本K-S检验或者正态检验时，有时会有错误提示“Kolmogorov - Smirnov检验里不应该有连结”，
#这是因为K-S检验只对连续CDF有效，而连续CDF中出现相同值的概率为0，也就是说数据中倘出现相同值，则连续分布的假设不成立，因此R会报错。这也提醒我们，在做正态性检验之前，要先对数据进行描述性分析，对数据整体要先有个大致的认识，这也才后续才能选择正确的检验方法。
#检验两个未知的分布是否为同一分布
X
2.36, 1.31, 1.05, -0.32, -0.40, 1.06, -2.47,0.39, 1.09, -1.28)
Y
ks.test(X,Y)
#p=0.5286>0.05,不能拒绝原假设，两个服从相同的分布
#注意：如果样本太少，可能出现错误结论
a1=c(-1020,-1040);a2=c(1020,1040)
t.test(a1,a2)
#一般情况，样本越小，往往容易得到不能拒绝原假设的结论
#pearson相关系数
a=c(1,3,5,7,9);b=c(1,4,6,9,10)
cor(a,b)
cor.test(a,b) #检验相关系数的显著性
cor(iris[1:4]) #相关系数,参数填数据集，则计算相关系数矩阵
a1=rnorm(5);b1=rnorm(5);cor(a1,b1);cor.test(a1,b1) #自己模拟生成两个变量
#spearman相关系数,亦即秩相关系数
#spearman和kendall都是等级相关系数，亦即其值与两个相关变量的具体值无关，而仅仅与其值之间的大小关系有关。
#spearman相关系数，亦即秩相关系数，根据随机变量的等级而不是其原始值衡量相关性的一种方法。
m=c(1,2,4,3);n=c(100,101,102,103)
m1=c(30,31,35,34);n1=c(85,87,90,93)
cor(m,n);cor(m1,n1)
cor(m,n,method = “spearman”);cor(m1,n1,method = “spearman”)
cor.test(m,n,method = “spearman”);cor.test(m1,n1,method = “spearman”)
#spearman相关系数的计算可以由计算pearson系数的方法,只需要把原随机变量中的原始数据替换成其在随机变量中的等级顺序即可:
acf #自相关和协方差函数
acf(airmiles,type=‘correlation’,lag.max=10) #自相关
pacf(airmiles,lag.max=10) #偏自相关
pairs(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=iris,
main=“Simple Scatterplot Matrix”) #散点图矩阵
install.packages(“scatterplot3d”) #3D散点图
library(scatterplot3d)
scatterplot3d(irisSepal.Length,irisSepal.Length, irisSepal.Length,irisPetal.Length, iris$Petal.Width)
install.packages(“corrgram”) #有兴趣的同学自己练习
library(corrgram)
#1、设置排序处理
corrgram(mtcars,order=TRUE)
#2、设置上下三角面板形状
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie)
#3、只显示下三角部分
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=NULL)
#4、调整面板颜色
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie,
col.regions=colorRampPalette(c(“darkgoldenrod4”,“burlywood1”,“white”, “darkkhaki”,“darkgreen”)))
install.packages(“corrplot”)
library(corrplot)
#1、使用不同的method绘制相关矩阵图
methods
par(mfrow=c(2,3))
t0=mapply(function(x){corrplot(cor(mtcars), method=x,order=“AOE”)},methods)
par(mfrow=c(1,1))
#2、设置method=color绘制热力矩阵图
corrplot(cor(mtcars), method=“color”, order = “AOE”,tl.col=“black”,tl.srt=45,
addCoef.col=“black”,col=colorRampPalette(c("#7F0000",“red”,"#FF7F00",
“yellow”,“white”, “cyan”, “#007FFF”, “blue”,"#00007F"))(20))
#3、绘制上下三角及不同色彩的相关矩阵图
library(RColorBrewer)
par(mfrow=c(2,2))
corrplot(cor(mtcars),type=“lower”)
corrplot(cor(mtcars),type=“lower”,order=“hclust”,
col=brewer.pal(n=8,name=“RdYlBu”))
corrplot(cor(mtcars),type=“upper”,order=“AOE”,
col=c(“black”,“white”),bg=“lightblue”)
corrplot(cor(mtcars),type=“upper”,order=“FPC”,
col=brewer.pal(n=8, name=“PuOr”))
par(mfrow=c(1,1))
d
hc
plot(hc)
rect.hclust(hc,k=3)
install.packages(“pvclust”)
library(pvclust)
cluster.bootstrap
plot(cluster.bootstrap)
pvrect(cluster.bootstrap) #自己练习部分结束
#回归分析
#回归分析的形式：fitted.model=lm(formula,data=data.frame)
#yx;yx+1 表示y对x的一元线性回归，默认包含常数项
y1=rnorm(30);x1=1:30
y2=3*x2-5;
x2=rnorm(30);k2=lm(y2~x2);summary(k2)
k1=lm(y1~x1)
summary(k1)
#log(y)~x1+x2 表示y的对数对自变量的回归，默认带常数项
#yploy(x,2),y1+x+I(x^2) 对x的一元二次多项式回归
k3=lm(y1~1+x1+I(x1^2));k3;summary(k3)
#y~A， 表示单因素方差分析，A是分类变量
#y~A+x 单因素协方差分析，A分类变量，x连续型变量
install.packages(“readxl”) #读取Excel数据的包
library(readxl)
y=c(16.3,16.8,19.2,18, 19.5,20.9,21.1,20.9,20.3,22,23.7,34.2,36.2)
x1=c(1.1,1.4,1.7,1.7,1.8,1.8,1.9,2,2.3,2.4,2.5,2.6,2.7)
x2=c(1.1,1.5,1.8,1.7,1.9,1.8,1.8,2.1,2.4,2.5,2.7,2.9,2.94)
A=c(0,1,1,0,0,1,1,1,0,0,2,2,2)
A=as.factor(A) #注意R语言中分类变量的处理
B=c(“wuhan”,“changsha”,“wuhan”,“hefei”,“changsha”,“wuhan”,“hefei”,“changsha”,“wuhan”,“hefei”,“wuhan”,“changsha”,“hefei”)
m1=lm(y~x1+x2+A)
m1
m2=aov(y~A) #单因素方差分析
m2
summary(m2)
m3=lm(y~x1+I(x1**2)) #多项式要用I括起来表示整体
m4=lm(y~x1+B)
m4
m5=lm(y~x1+B+x1:B);m5 #交互项用：
m6=lm(y~.,data=) #数据集中除因变量的所有其它变量作为自变量
m7=lm(y~x1x2B);m7 #三变量单独、两两交互、三者交互的组合
m8=lm(y~(x1+x2+B)**2);m8 #单独、两两交互
#summary(详细结果)coefficients(模型参数)confint(置信区间)
#fitted(拟合值)residuals(残差值)AIC(赤池信息)predict(预测值)
m9=lm(y~x1+x2)
summary(m9)
coefficients(m9)
confint(m9)
AIC(m9)
fitted(m9)
residuals(m9)
predict(m9)
formula(m3)
plot(x1,y) #散点图
boxplot(x1~B) #箱线图
plot(m3) #默认画4个图，残差对拟合值图，残差正态概率图，标准化残差绝对值平方根对拟合值图，cook距离图
k1=par(mfrow=c(1,1)) #图的展示方式，2行2列
plot(m3)
k2=lm(y~x1)
plot(y,x1)
abline(lm(y~x1)) #散点图加拟合线
#回归模型诊断
install.packages(“car”) #很多回归诊断功能
library(car)
library(carData)
durbinWatsonTest(m3) #自相关检验
crPlots(m3) #看线性趋势
ncvTest(m3) #原假设是同方差
vif(m3) #方差膨胀因子
outlierTest(m3) #原假设无离群点