corpcor程序包 自行下载加载(通过界面上的 程序包)
复制下面文本到R语言控制台(console)中 粘贴 回车 可以看到示例结果
以下内容转自普蘭塔论坛 作者gagea
- #数量化方法I
- library(corpcor)
- y=c(3,5,6,7,9,11,9,7,7,6)
- c11=c(1,0,0,1,0,0,1,0,0,0)
- c12=c(0,1,0,0,1,0,0,1,0,1)
- c13=c(0,0,1,0,0,1,0,0,1,0)
- c21=c(1,1,1,0,0,0,0,0,1,1)
- c22=c(0,0,0,1,1,1,1,1,0,0)
- test=data.frame(y=y,c11=c11,c12=c12,c13=c13,c21=c21,c22=c22)
- #定义数量化理论I分析函数
- ##x为反应矩阵,用数据框表示;y为基准变量,用向量表示;nItem为项目数;nCati为各项目类目数,用向量表示。
- ##ItemNames为项目名称向量。
- TOQO=function(x,y,nItem,ItemNames,nCati){
- xnames=names(x)
- xm=as.matrix(x)
- ym=as.matrix(y)
- x.x=t(xm)%*%xm#正则方程系数矩阵
- x.y=t(xm)%*%ym#正则方程右端常数向量
- #对每个j=2,...,m,删去第j项目第一类目的方程
- Index=1:sum(nCati)
- IndexDel=NULL
- for(i in 2:nItem) {Index=Index[Index!=sum(nCati[1:(i-1)])+1]#确定保留索引号
- IndexDel[i-1]=sum(nCati[1:(i-1)])+1#确定要删除的索引号
- }
- x.x_2=x.x[Index,Index]
- x.y_2=x.y[Index]
- b1=solve(x.x_2,x.y_2)#预测方程部分解
- b=NULL#初始化预测方程完整解。
- b[Index]=b1
- b[IndexDel]=0
- ItemM=matrix(nrow=nrow(x),ncol=nItem)#项目预测得分矩阵
- for(i in 1:nItem){
- if(i==1) ItemM[,i]=xm[,1:nCati[i]]%*%b[1:nCati[i]]
- else ItemM[,i]=xm[,IndexDel[i-1]:(IndexDel[i-1]+nCati[i]-1)]%*%b[IndexDel[i-1]:(IndexDel[i-1]+nCati[i]-1)]
- }
- M=cbind(ym,ItemM)
- attributes(M)$dimnames[[2]]=c("y",ItemNames)
- xcor=cor(M)#相关系数矩阵
- xpcor=cor2pcor(xcor)#偏相关矩阵
- pcor=xpcor[2:(nItem+1),1]#偏相关系数
- pcordf=as.data.frame(t(pcor))
- ItemVar=apply(ItemM,2,var)
- yVar=var(ym)
- ItemVarRatio=ItemVar/yVar#方差比
- IVR=as.data.frame(t(ItemVarRatio))
- Range=NULL
- for(i in 1:nItem){#范围
- if(i==1) Range[i]=max(b[1:nCati[i]])-min(b[1:nCati[i]])
- else Range[i]=max(b[IndexDel[i-1]:(IndexDel[i-1]+nCati[i]-1)])-min(b[IndexDel[i-1]:(IndexDel[i-1]+nCati[i]-1)])
- }
- Range=as.data.frame(t(Range))
- Table=rbind(pcordf,IVR,Range)
- names(Table)=ItemNames
- row.names(Table)=c("偏相关系数","方差比","范围")
- ans=list(b=b,Table=Table)
- ans
- }
- #例:
- x=test[,2:6]
- y=test[,1]
- nItem=2
- nCati=c(3,2)
- ItemNames=c("体重","性别")
- TOQO(x,y,nItem,ItemNames,nCati)
