利用R求简单相关系数及其显著性检验
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
题目:
测定“丰产三号”小麦的每株穗数x1,主茎上每穗结实小穗数x2,百粒重x3(单位:g),主茎株高x4(单位:cm)和每株籽粒产量y(单位:g)的15
丰产3号小麦栽培试验的观测值
利用R软件自带函数计算简单相关系数及其检验p值:
手工编写相关系数计算函数和检验函数:
COR.test = function(X,R); #求F检验的p值,为矩阵形式
COR = function(X); #求相关系数,为矩阵形式
lxy = function(x,y); #各个观测值之间的离均差乘积
得到的简单相关系数表:
与R软件自带的函数得到的结果相比较,可以看出相关系数结果一致,检验p 值偏差可以忽略。
各个函数具体代码如下:
COR.test = function(X,R){options(digits = 4) #求F检验的p值,为矩阵形式n = dim(X)[2]; #得到p矩阵的阶数
p = diag(0,n) ; #n阶零矩阵
for (i in 1:n)
{
for (j in 1:n){
f = R[i,j]^2/((1-R[i,j]^2)/(dim(X)[1]-2)); #用F检验对相关系数作显著性检验
p[i,j]=1-pf(f,1,dim(X)[1]-2); #用F检验计算p值
}
}
p
}
COR = function(X){options(digits = 3) #求相关系数,为矩阵形式
n = dim(X)[2]; #得到相关系数矩阵的阶数
Y = diag(n); #产生单位矩阵
for (i in 1:n)
for( j in 1:n)
Y[i,j] = lxy(X[,i],X[,j])/sqrt(lxy(X[,i],X[,i])*lxy(X[,j],X[,j]));#简单相关系数计算公式
Y
}
lxy = function(x,y){ #各个观测值之间的离均差乘积
n = length(x)
sum(x*y) - sum(x)*sum(y)/n
}
f_COR = function(R){ #复相关系数,R为简单相关系数矩阵R = solve(R) #简单相关系数矩阵的逆矩阵
n = dim(R)[1]
r = numeric(n) #产生n阶零向量
for (i in 1:n){
r[i] = sqrt(1 - 1/R[i,i]) 复相关系数的计算公式
}
r
}
p_COR = function(X,R){ #偏相关系数,X为输入数据框或者矩阵,R为简单相关系数矩阵
R = solve(R) #简单相关系数矩阵的逆矩阵
n = dim(R)[1]
n1 = dim(X)[1]
p = dim(X)[2] - 1 #自由度为n1-p-1
r = diag(0,n) #产生n阶零矩阵
for (i in 1:(n - 1)){ #对角线以上部分
for (j in (i + 1):n){
r[i,j] = -(1)*R[i,j]/sqrt(R[i,i]*R[j,j]) #偏相关系数计算公式}
}
for (i in 1:(n - 1)){
for (j in (i + 1):n){ #对角线以下部分
F = r[i,j]^2/(1 - r[i,j]^2)*( n1 - p - 1 ) #F检验值
r[j,i] =1 - pf(F,1,n1 - p - 1,ncp = 0,log = FALSE) #F检验得到的p值}
}
r #对角线以上为偏相关系数,以下为其检验的p值
}