如何在Matlab中复现R语言chisq.test(cbind(x,y))的卡方检验结果

我在R和Matlab中对两个观测频率数据集x、y做列联表卡方检验，尝试复现一致结果，重点关注x=y的特殊场景：

• 在R中，执行chisq.test(x,y)得到卡方值12、自由度9、p值0.2133；但执行chisq.test(cbind(x,y))时，因x=y，得到符合直觉的卡方值0、自由度3、p值1。
• 在Matlab中，[~,chi2,p] = crosstab(x,y)的结果和R里的chisq.test(x,y)一致，但无法复现chisq.test(cbind(x,y))的结果。我自行编写代码复现crosstab逻辑，仍未找到实现目标的方法。

如何在Matlab中实现当x=y时，卡方检验统计量为0、p值为1的结果，即复现R语言chisq.test(cbind(x,y))的功能？

我编写了如下Matlab代码复现crosstab(x,y)的逻辑，但未实现目标功能：

% Inputs
x = [1000 100 10 1];
y = [1000 100 10 1];
% Chi-squared test on contingency tables
[unique_x,ix,jx] = unique(x);
[unique_y,iy,jy] = unique(y);
OT = zeros(numel(unique_x),numel(unique_y));
for i = 1:numel(x)
     OT(jx(i),jy(i)) = OT(jx(i),jy(i)) + 1;
end
R = sum(OT,2);
C = sum(OT);
n = sum(sum(OT));
ET = (R*C)./n;
chi2 = (OT - ET).^ 2 ./ ET;
chi2 = sum(chi2(:))                  % <-- test statistic
df = (numel(R)-1)*(numel(C)-1)       % <-- degrees of freedom
p = gammainc(chi2/2,df/2,'upper')    % <-- p-value
OT                                   % <-- OT = Observed frequency contingency table 
ET                                   % <-- ET = Expected frequency contingency table

代码返回结果：

chi2 =
    12
df =
     9
p =
      0.21331
OT =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
ET =
         0.25         0.25         0.25         0.25
         0.25         0.25         0.25         0.25
         0.25         0.25         0.25         0.25
         0.25         0.25         0.25         0.25

要复现R中chisq.test(cbind(x,y))的功能，本质是把x和y直接作为列联表的两列（每行对应同一类别在两个变量中的频数），执行多分类配对卡方检验，而非将x、y作为两个分类变量构建交叉表。

在Matlab中可直接构造目标列联表并计算，代码如下：

% 输入数据
x = [1000 100 10 1];
y = [1000 100 10 1];

% 构造列联表：每行对应一个类别，两列分别为x、y的观测频数
contingency_table = [x', y'];

% 计算期望频数：假设类别在x、y中无差异，每行两列的期望为该行总和的一半
row_sums = sum(contingency_table, 2);
ET = [row_sums/2, row_sums/2];

% 计算卡方统计量
chi2 = sum(sum((contingency_table - ET).^2 ./ ET));

% 自由度：类别数 - 1（多分类配对卡方的自由度为类别数减1）
df = size(contingency_table, 1) - 1;

% 计算p值
p = gammainc(chi2/2, df/2, 'upper');

% 输出结果
disp(['卡方值: ', num2str(chi2)]);
disp(['自由度: ', num2str(df)]);
disp(['p值: ', num2str(p)]);

运行后会得到和R中一致的结果：

卡方值: 0
自由度: 3
p值: 1

关键说明 #

• 原crosstab逻辑是统计x、y各分类组合的出现次数，对应R的chisq.test(x,y)；而我们需要的是将x、y作为列联表的两列，直接对比同一类别在两个变量中的频数。
• 当x=y时，观测矩阵与期望矩阵完全一致，卡方值自然为0，p值为1。

内容的提问来源于stack exchange，提问作者BostonPlummer