定性数据分析中高维列联表分析流程（附R代码）

第一篇

关于定性数据分析

上学期选修了定性数据分析这门课程，教材是王静龙、梁小筠的定性数据分析。

高维列联表分析流程

目录
第一章 绪论 7
1.1 问题研究背景 7
1.2 数据来源 7
1.3 研究意义 1
第二章 高维列联表独立性检验 2
2.1 高维列联表的相互独立性检验 2
2.2 高维列联表的边缘独立性检验 2
2.3 高维列联表的条件独立性检验 3
2.3.1层属性“地区”给定时的条件独立性检验 3
2.3.2行属性“性别”给定时的条件独立性检验 4
第三章 高维列联表相合性检验 5
3.1 高维列联表的压缩相合性检验 5
3.2 高维列联表分层相合性检验 5
3.3 高维列联表条件相合性检验（Cochran-Mantel-Haenszel检验） 7
3.2 高维列联表相合齐次性检验（Breslow-Day检验） 8
第四章 结论与建议 10
4.1 结论 10
4.2 相关建议 11
参考文献 12

直接上代码给需要的人（R）

// An highlighted block
附录4
###################定性数据统计分析论文##########
#####三维列联表的压缩与分层
##三维列联表压缩函数，输入为三维列联表(行*列*层)数组，按层压缩，得到边缘表
compress=function(x){
  x_compress=matrix(0,nrow = dim(x)[1],ncol = dim(x)[2])#构建一个与边缘表相同大小的0矩阵
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    x_compress=x_compress+x[,,k]#x[,,k]第k层数据矩阵
  }##dim(z)[3]#数组z的层数
  return(x_compress)#返回边缘表
}

##三维列联表分层与相合性检验函数，输入为三维列联表(行*列*层)数组，按层分层，得到部分表和检验结果
layering=function(x){
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    x_layer=x[,,k]#x[,,k]第k层数据矩阵
    print(x_layer)#打印部分表
    congruence_test(x_layer,alternative="greater")###相合性检验，H1:正相合
    cat('\n')
  }##dim(z)[3]#数组z的层数
}
congruence_test=function(x,alternative="")
#适用于列联表的相合性度量与检验问题
#x为列联表矩阵,alternative对应于备择假设
{
  n=sum(x)
  G=0;H=0
  r=nrow(x)
  c=ncol(x)
  r1=r-1;c1=c-1
  #计算G:
  for (i in 1:r1){
    for (j in 1:c1){
      G=G+x[i,j]*sum(x[(i+1):r,(j+1):c])
    }
  }
  #计算H：
  for (i in 1:r1){
    for (j in 2:c){
      H=H+x[i,j]*sum(x[(i+1):r,1:(j-1)])
    }
  }
  #计算z,TA,TB,TAB:
  z=G-H
  TA=sum(rowSums(x)*(rowSums(x)-1)/2)
  TB=sum(colSums(x)*(colSums(x)-1)/2)
  #TAB=G+H+TA+TB-n*(n-1)/2
  Cn2=n*(n-1)/2
  ######相合性度量######
  #计算各系数的值
  Kendall_tau=z/sqrt((Cn2-TA)*(Cn2-TB))###Kendall-tau系数
  Gamma=(G-H)/(G+H)######Gamma系数
  #####Somers-d系数：
  d_BA=(G-H)/(Cn2-TA)###列联表为2*c,列属性B依赖于行属性A，用d_B|A度量相合性
  d_AB=(G-H)/(Cn2-TB)###列联表为r*2,行属性A依赖于列属性B，用d_A|B度量相合性
  
  #######相合性检验######
  #sigma平方的近似计算公式,表示sigma的平方
  sigma_2=(n^3-sum(rowSums(x)^3))*(n^3-sum(colSums(x)^3))/(9*n^3) 
  #构建检验统计量：U统计量
  U=z/sqrt(sigma_2)
  if(alternative=="two.side")
  {p_value=1-pchisq(U^2, 1)}
  else 
  {
    if(alternative=="greater")####greater：备择假设为正相合
    {p_value=pnorm(-U)}
    else if(alternative=="less")####less:备择假设为负相合
    {p_value=pnorm(U)}
    else{cat("please input:\n alternative= 'two.side','greater',or'less'")}
  }
  cat('关于相合性度量的相关系数：\n')
  cat('G=',G,'\n')
  cat('H=',H,'\n')
  cat('z=',z,'\n')
  cat('T_A=',TA,'\n')
  cat('T_B=',TB,'\n')
  cat('Kendall_tau=',Kendall_tau,'\n')
  cat('Gamma=',Gamma,'\n')
  cat('d_B|A=',d_BA,'\n')
  cat('d_A|B=',d_AB,'\n\n')
  cat('相合性检验:\n')
  cat('检验统计量的值U=',U,'\n')
  cat('p_value=',p_value,'\n')
}
##高维(三维)列联表的条件独立性检验
Conditional_independence_test=function(x){
  Likelihood_ratio_layer=0
  Likelihood_ratio_row=0
  Likelihood_ratio_col=0
  #层属性C给定时计算似然比检验统计量的值
  #原假设H0：层属性C给定后行属性A与列属性B条件独立，备择假设H1：层属性C给定后行属性A与列属性B不条件独立
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    for (i in 1:dim(x)[1]) {##dim(z)[3]#数组z的层数
      for (j in 1:dim(x)[2]) {
        Likelihood_ratio_layer=Likelihood_ratio_layer+(-2*x[i,j,k]*log(sum(x[i,,k])*sum(x[,j,k])/(sum(x[,,k])*x[i,j,k])))
      }
    }
  }
  #层属性C给定时计算条件独立性检验的p值
  p_value_layer=1-pchisq(Likelihood_ratio_layer, dim(x)[3]*(dim(x)[1]-1)*(dim(x)[2]-1))
  
  #行属性A给定时计算似然比检验统计量的值
  #原假设H0：行属性A给定后层属性C与列属性B条件独立，备择假设H1：行属性A给定后层属性C与列属性B不条件独立
  for (i in 1:dim(x)[1]) {
    for (j in 1:dim(x)[2]) {
      for (k in 1:dim(x)[3]) {
        Likelihood_ratio_row=Likelihood_ratio_row+(-2*x[i,j,k]*log(sum(x[i,,k])*sum(x[i,j,])/(sum(x[i,,])*x[i,j,k])))
      }
    }
  }
  #行属性A给定时计算条件独立性检验的p值
  p_value_row=1-pchisq(Likelihood_ratio_row, dim(x)[1]*(dim(x)[2]-1)*(dim(x)[3]-1))
  
  #列属性B给定时计算似然比检验统计量的值
  #原假设H0:列属性B给定后行属性A与层属性C条件独立，备择假设H1：列属性B给定后行属性A与层属性C不条件独立
  for (j in 1:dim(x)[2]) {
    for (i in 1:dim(x)[1]) {
      for (k in 1:dim(x)[3]) {
        Likelihood_ratio_col=Likelihood_ratio_col+(-2*x[i,j,k]*log(sum(x[i,j,])*sum(x[,j,k])/(sum(x[,j,])*x[i,j,k])))
      }
    }
  }
  #列属性B给定时计算条件独立性检验的p值
  p_value_col=1-pchisq(Likelihood_ratio_col, dim(x)[2]*(dim(x)[1]-1)*(dim(x)[3]-1))
  
  cat('层属性给定时条件独立性检验的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Layer=',Likelihood_ratio_layer,'\n')
  cat('P_value-Layer=',p_value_layer,'\n')
  cat('行属性给定时条件独立性检验的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Row=',Likelihood_ratio_row,'\n')
  cat('P_value-Row=',p_value_row,'\n')
  cat('列属性给定时条件独立性检验的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Col=',Likelihood_ratio_col,'\n')
  cat('P_value-Col=',p_value_col,'\n')
  
}
##高维列联表的独立性检验情况一(相互独立-H0:属性A、B、C相互独立<->H1:属性A、B、C不相互独立)
Independent_test=function(x){
  Likelihood_ratio_eachother=0
  #原假设H0：属性A、B、C相互独立，备择假设H1：属性A、B、C不相互独立
  #计算似然比检验统计量
  for (i in 1:dim(x)[1]) {
    for (j in 1:dim(x)[2]) {
      for (k in 1:dim(x)[3]) {
        Likelihood_ratio_eachother=Likelihood_ratio_eachother+(-2*x[i,j,k]*log(sum(x[i,,])*sum(x[,j,])*sum(x[,,k])/((sum(x[,,]))^2*x[i,j,k])))
      }
    }
  }
  #计算属性A、B、C的相互独立的独立性检验的p值
  p_value_eachother=1-pchisq(Likelihood_ratio_eachother, dim(x)[1]*dim(x)[2]*dim(x)[3]-dim(x)[1]-dim(x)[2]-dim(x)[3]+2)
  
  cat('属性A、B、C相互独立的独立性检验的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Eachother=',Likelihood_ratio_eachother,'\n')
  cat('P_value-Eachother=',p_value_eachother,'\n')
}

##高维列联表的独立性检验情况二(边缘独立性检验)
#原假设H0:属性A和(B,C)相互独立<->备择假设H1:A和(B,C)不相互独立,属性B和(A,C)及属性C和(A,B)类似
Edge_independence_test=function(x){
  Likelihood_ratio_layer_edge=0
  Likelihood_ratio_row_edge=0
  Likelihood_ratio_col_edge=0
  #层属性C与另外俩个属性合并所得的边缘表
  #计算似然比检验统计量的值
  #原假设H0：层属性C和(A,B)相互独立<->备择假设H1:C和(A,B)不相互独立
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    for (i in 1:dim(x)[1]) {##dim(z)[3]#数组z的层数
      for (j in 1:dim(x)[2]) {
        Likelihood_ratio_layer_edge=Likelihood_ratio_layer_edge+(-2*x[i,j,k]*log(sum(x[,,k])*sum(x[i,j,])/(sum(x[,,])*x[i,j,k])))
      }
    }
  }
  #层属性C与另外俩个属性合并所得的边缘表下计算边缘独立性检验的p值
  p_value_layer_edge=1-pchisq(Likelihood_ratio_layer_edge, (dim(x)[3]-1)*(dim(x)[1]*dim(x)[2]-1))
  
  #行属性A与另外俩个属性BC合并所得的边缘表
  #计算似然比检验统计量的值
  #原假设H0：层属性A和(B,C)相互独立<->备择假设H1:A和(B,C)不相互独立
  for (i in 1:dim(x)[1]) {
    for (j in 1:dim(x)[2]) {
      for (k in 1:dim(x)[3]) {
        Likelihood_ratio_row_edge=Likelihood_ratio_row_edge+(-2*x[i,j,k]*log(sum(x[i,,])*sum(x[,j,k])/(sum(x[,,])*x[i,j,k])))
      }
    }
  }
  #行属性A与另外俩个属性BC合并所得的边缘表的边缘独立性检验的p值
  p_value_row_edge=1-pchisq(Likelihood_ratio_row_edge, (dim(x)[1]-1)*(dim(x)[2]*dim(x)[3]-1))
  
  #列属性B与另外俩个属性AC合并所得的边缘表
  #计算似然比检验统计量的值
  #原假设H0：列属性B和(A,C)相互独立<->备择假设H1:B和(A,C)不相互独立
  for (j in 1:dim(x)[2]) {
    for (i in 1:dim(x)[1]) {
      for (k in 1:dim(x)[3]) {
        Likelihood_ratio_col_edge=Likelihood_ratio_col_edge+(-2*x[i,j,k]*log(sum(x[,j,])*sum(x[i,,k])/(sum(x[,,])*x[i,j,k])))
      }
    }
  }
  #列属性B与另外俩个属性AC合并所得的边缘表的边缘独立性检验的p值
  p_value_col_edge=1-pchisq(Likelihood_ratio_col_edge, (dim(x)[2]-1)*(dim(x)[1]*dim(x)[3]-1))
  
  cat('层属性C与另外俩个属性AB合并下所得的边缘表的相互独立性检验(边缘表独立性检验)的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Layer_Edge=',Likelihood_ratio_layer_edge,'\n')
  cat('P_value-Layer_Edge=',p_value_layer_edge,'\n')
  cat('行属性A与另外俩个属性BC合并下所得的边缘表的相互独立性检验(边缘表独立性检验)的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Row_Edge=',Likelihood_ratio_row_edge,'\n')
  cat('P_value-Row_Edge=',p_value_row_edge,'\n')
  cat('列属性B与另外俩个属性AC合并下所得的边缘表的相互独立性检验(边缘表独立性检验)的似然比检验统计量及p值为：\n')
  cat('Likelihood_Ratio_Col_Edge=',Likelihood_ratio_col_edge,'\n')
  cat('P_value-Col_Edge=',p_value_col_edge,'\n')
}

##条件相合性检验(Cochran-Mantel-Haenszel chi-square检验)
#H0：层属性C给定后行属性A与列属性B条件独立<->H1:层属性C给定后行属性A与列属性B条件正向合(条件负向合、条件相合)
Cochran_Mantel_Haenszel=function(x,alternative=""){#条件相合性检验
  U_test_statistic=0
  U_part_up1=0
  U_part_up2=0
  U_part_low=0
  
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    U_part_up1=U_part_up1+x[1,1,k]#公式上部分前部分
    U_part_up2=U_part_up2+sum(x[1,,k])*sum(x[,1,k])/sum(x[,,k])#公式上部分后部分
    U_part_low=U_part_low+sum(x[1,,k])*sum(x[2,,k])*sum(x[,1,k])*sum(x[,2,k])/((sum(x[,,k]))^2*(sum(x[,,k])-1))#公式下部分未开方部分
  }
  U_test_statistic=(U_part_up1-U_part_up2)/sqrt(U_part_low)
  chi_square_statistic=U_test_statistic^2
  
  if(alternative=="two.side")####two.side：备择假设为条件相合
  {p_value=1-pchisq(chi_square_statistic, 1)}
  else 
  {
    if(alternative=="greater")####greater：备择假设为条件正相合
    {p_value=pnorm(-U_test_statistic)}
    else if(alternative=="less")####less:备择假设为条件负相合
    {p_value=pnorm(U_test_statistic)}
    else{cat("please input:\n alternative= 'two.side','greater',or'less'")}
  }
  cat('条件相合性检验:\n')
  cat('检验统计量的值U_test_statistic=',U_test_statistic,'\n')
  cat('p_value=',p_value,'\n')
}


###高维列联表条件相合齐性检验(Breslow-Day chi-square检验)(优比检验)
#检验各层条件相合程度是否相同
#原假设H0:各层条件相合程度相同theta_1=...=theta_t<->H1:各层相合程度不全相同theta_1...theta_t不全相等
Breslow_Day_test=function(x){
  theta_estimate=0
  eta_mean_up=0
  eta_mean_low=0
  Chi_square_statistic=0
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    theta_estimate=x[1,1,k]*x[2,2,k]/(x[1,2,k]*x[2,1,k])
    a=1/x[1,1,k]+1/x[1,2,k]+1/x[2,1,k]+1/x[2,2,k]
    eta_k=log(theta_estimate)
    eta_mean_up=eta_mean_up+eta_k/a
    eta_mean_low=eta_mean_low+1/a
  }
  eta_mean=eta_mean_up/eta_mean_low
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    theta_estimate=x[1,1,k]*x[2,2,k]/(x[1,2,k]*x[2,1,k])
    a=1/x[1,1,k]+1/x[1,2,k]+1/x[2,1,k]+1/x[2,2,k]
    eta_k=log(theta_estimate)
    Chi_square_statistic=Chi_square_statistic+(eta_k-eta_mean)^2/a
  }
  p_value=1-pchisq(Chi_square_statistic,dim(x)[3]-1)
  cat('各层条件相合齐性检验统计量及p值为：\n')
  cat('Chi_Square_Statistic=',Chi_square_statistic,'\n')
  cat('P_value=',p_value,'\n')
}


##导入第四章相合性检验函数，用于计算边缘表和部分表的相合性检验
source("D://Rcode//定性数据统计分析//第四章//第4章.r",encoding = "utf-8")

##导入第五章
source("D://Rcode//定性数据统计分析//第五章//第5章.r",encoding = "utf-8")

data=read.csv('三维列联表csv格式.csv')
print(data)
data1=data[-1,-1]
print(data1)
data2=t(data1)
data3=as.vector(data2)
print(data3)
data3=data3/10000

dim1 <- c("男","女")
dim2 <- c("未上学","小学","初中","高中","专科","本科","研究生")#列属性
dim3 <- c("北京","天津","河北","山西","内蒙古","辽宁","吉林","黑龙江","上海","江苏","浙江","安徽","福建","江西","山东","河南","湖北","湖南","广东","广西","海南","重庆","四川","贵州","云南","西藏","陕西","甘肃","青海","宁夏","新疆")#层属性
##myarray <- array(vector, dimensions, dimnames)
population <- array(data3,c(2,7,31),dimnames = list(dim1,dim2,dim3))#三维4*4*3列联表
population

##独立性检验
Independent_test(population)

##边缘独立性检验
Edge_independence_test(population)

##高维列联表的独立性检验情况三(条件独立性检验)
Conditional_independence_test(population)
################三维列联表的压缩与检验
z_compress=compress(population)
z_compress#压缩后所得边缘表
#对压缩后所得边缘表进行相合性检验
congruence_test(z_compress,alternative="less")

layering1=function(x){
  for (k in 1:dim(x)[3]) {##dim(z)[3]#数组z的层数
    x_layer=x[,,k]#x[,,k]第k层数据矩阵
    print(k)
    print(x_layer)#打印部分表
    
    congruence_test(x_layer,alternative="less")###相合性检验，H1:负相合
    cat('\n')
  }##dim(z)[3]#数组z的层数
}

#分层相合性检验
layering1(population)

#条件相合性检验
Cochran_Mantel_Haenszel(population,alternative = "less")

#相合其次性检验
#Breslow_Day_test(population)

dim4 <- c("男","女")
dim5 <- c("未上学","上学")#列属性
dim6 <- c("北京","天津","河北","山西","内蒙古","辽宁","吉林","黑龙江","上海","江苏","浙江","安徽","福建","江西","山东","河南","湖北","湖南","广东","广西","海南","重庆","四川","贵州","云南","西藏","陕西","甘肃","青海","宁夏","新疆")#层属性

for (k in 1:dim(population)[3]) {##dim(z)[3]#数组z的层数
  population1=population[,,k]#x[,,k]第k层数据矩阵
  print(k)
  a=population1[1,1]
  b=population1[1,2]+population1[1,3]+population1[1,4]+population1[1,5]+population1[1,6]+population1[1,7]
  c=population1[2,1]
  d=population1[2,2]+population1[2,3]+population1[2,4]+population1[2,5]+population1[2,6]+population1[2,7]
  population2<- array(c(a,c,b,d),c(2,2,31),dimnames = list(dim4,dim5,dim6))#
  print(population2)#打印部分表
  
  cat('\n')
}##dim(z)[3]#数组z的层数

Breslow_Day_test(population2)

关注微信公众号：菜田里的守望者
回复：定性数据分析 获取更多定性数据分析相关源码

需要完整代码或书籍PDF版的小伙伴可关注微信公众号：菜田里守望者

打开微信扫一扫关注吧，你们的支持就是我的动力