R的典型相关分析感觉比较复杂,还是SAS简单

数据来自SAS:



Weight    Waist    Pulse    Chins    Situps    Jumps

191    36    50    5    162    60

189    37    52    2    110    60

193    38    58    12    101    101

162    35    62    12    105    37

189    35    46    13    155    58

182    36    56    4    101    42

211    38    56    8    101    38

167    34    60    6    125    40

176    31    74    15    200    40

154    33    56    17    251    250

169    34    50    17    120    38

166    33    52    13    210    115

154    34    64    14    215    105

247    46    50    1    50    50

193    36    46    6    70    31

202    37    62    12    210    120

176    37    54    4    60    25

157    32    52    11    230    80

156    33    54    15    225    73

138    33    68    2    110    43


<br />
fit <- read.table("D:/data/fit.txt", header=TRUE, sep="\t", na.strings="0.0", dec=".", strip.white=TRUE)<br />
attach(fit)<br />
#空间统计包<br />
library(fields)<br />
# t() 转置函数<br />
#stats()  fields包中的函数,显示数据的一些描述性统计量<br />
t(stats(fit))<br />
<br />
# 将变量分成两类<br />
wwp<-fit[,1:3]<br />
csj<-fit[,4:6]<br />
#cancor()在stats包中,典型相关分析<br />
cancor(wwp,csj)<br />
#CCA 典型相关分析包<br />
library(CCA)<br />
# 计算两类数据的相关矩阵 ,输出结果也包括但数据内变量相关系数<br />
matcor(wwp,csj)<br />
# cc 增强典型相关分析<br />
cc1 <- cc(wwp,csj)<br />
# 显示相关性<br />
cc1[1]<br />
# raw canonical coefficients <br />
 cc1[3:4]<br />
# 进一步分析,comput为cc的补充<br />
 cc2<-comput(wwp,csj, cc1)<br />
#显示典型变量及其相关系数<br />
 cc2[3:6]<br />
 <br />
# 标准维数检验<br />
<br />
ev<-cc1$cor^2<br />
ev2<-1-ev<br />
n<-dim(wwp)[1]<br />
p<-length(wwp)<br />
q<-length(csj)<br />
m<-n -3/2 - (p+q)/2<br />
w<-cbind(NULL)  # initialize wilks lambda<br />
<br />
for (i in 1:3){<br />
    w<-cbind(w,prod(ev2[i:3]))<br />
  }<br />
<br />
d1<-cbind(NULL)<br />
d2<-cbind(NULL)<br />
f<-cbind(NULL)  # initialize f<br />
for (i in 1:3){<br />
    s<-sqrt((p^2*q^2-4)/(p^2+q^2-5))<br />
    si<-1/s<br />
    df1<-p*q<br />
    d1<-cbind(d1,df1)<br />
    df2<-m*s-p*q/2+1<br />
    d2<-cbind(d2,df2)<br />
    r<-(1-w[i]^si)/w[i]^si<br />
    f<-cbind(f,r*df2/df1)<br />
    p<-p-1<br />
    q<-q-1<br />
}<br />
pv<-pf(f,d1,d2,lower.tail=FALSE)    #F分布<br />
dmat<-cbind(t(w),t(f),t(d1),t(d2),t(pv))<br />
colnames(dmat)<-c("WilksL","F","df1","df2","p")<br />
rownames(dmat)<-c(seq(1:length(w)))<br />
dmat<br />
#  wwp标准化典型相关系数<br />
<br />
sd<-sd(wwp)<br />
s1<-diag(sd)       # diagonal matrix of wwp sd's<br />
s1 %*% cc1$xcoef     # %*% 矩阵乘法<br />
sd<-sd(csj)<br />
s2<-diag(sd)       # diagonal matrix of acad sd's<br />
s2 %*% cc1$ycoef<br />
<br />




 WilksL          F df1      df2          p1

0.3503905 2.04823353   9 34.22293 0.063530942

0.9547227 0.17578229   4 30.00000 0.949120253

0.9947336 0.08470926   1 16.00000 0.77475327