IMO 最简单的方法是自己构建设计矩阵，因为它glmfit接受原始（观察到的）值矩阵或设计矩阵。一旦你编写了完整的模型，编写一个交互项就不是那么困难了。假设我们有两个预测变量，（连续）和（分类，具有三个无序级别，例如）。使用 Wilkinson 的符号，我们将把这个模型写成，忽略左侧（对于二项式结果，我们将使用 logit 链接函数）。我们只需要两个虚拟向量来编码级别（特定观察的存在/不存在），所以我们将有 5 个回归系数，加上一个截距项。这可以概括为$x$$g$$g={1,2,3}$y ~ x + g + x:gg

$$\beta_0 + \beta_1\cdot x +\beta_2\cdot\mathbb{I}_{g=2} +\beta_3\cdot\mathbb{I}_{g=3} + \beta_4\cdot x\times\mathbb{I}_{g=2} + \beta_5\cdot x\times\mathbb{I}_{g=3},$$

其中代表编码水平的指示矩阵。$\mathbb{I}$$g$

在 Matlab 中，使用在线示例，我将执行以下操作：

x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';
g = [1 1 1 1 2 2 2 2 3 3 3 3]';
gcat = dummyvar(g);
gcat = gcat(:,2:3); % remove the first column
X = [x gcat x.*gcat(:,1) x.*gcat(:,2)]; 
n = [48 42 31 34 31 21 23 23 21 16 17 21]';
y = [1 2 0 3 8 8 14 17 19 15 17 21]';
[b, dev, stats] = glmfit(X, [y n], 'binomial', 'link', 'probit');

我没有为拦截包含一列，因为它默认包含在内。设计矩阵看起来像

    2100           0           0           0           0
    2300           0           0           0           0
    2500           0           0           0           0
    2700           0           0           0           0
    2900           1           0        2900           0
    3100           1           0        3100           0
    3300           1           0        3300           0
    3500           1           0        3500           0
    3700           0           1           0        3700
    3900           0           1           0        3900
    4100           0           1           0        4100
    4300           0           1           0        4300

并且您可以看到交互项只是编码为与（g=2 和 g=3，因为我们不需要第一级）x的对应列的乘积。g

结果如下所示，作为系数、标准误差、统计量和 p 值（来自stats结构）：

   int.   -3.8929    2.0251   -1.9223    0.0546
   x       0.0009    0.0008    1.0663    0.2863
   g2     -3.2125    2.7622   -1.1630    0.2448
   g3     -5.7745    7.5542   -0.7644    0.4446
   x:g2    0.0013    0.0010    1.3122    0.1894
   x:g3    0.0021    0.0021    0.9882    0.3230

现在，可以通过计算与上述完整模型和简化模型（省略交互项，即设计矩阵的最后两列）的偏差差异来测试交互。这可以手动完成，也可以使用lratiotest提供似然比假设检验的功能。完整模型的偏差是 4.3122 ( dev)，而没有交互作用的模型是 6.4200（我用过glmfit(X(:,1:3), [y n], 'binomial', 'link', 'probit');），并且相关的 LR 检验有两个自由度（两个模型之间参数数量的差异）。由于比例偏差只是 GLM 的对数似然的两倍，我们可以使用

[H, pValue, Ratio, CriticalValue] = lratiotest(4.3122/2, 6.4200/2, 2)

其中统计量分布为（临界值为 5.9915，请参阅）。输出表明一个不显着的结果：我们不能断定观察到的样本之间存在相互作用。$\chi^2$chi2inv(0.95, 2)xg

我想你可以用你选择的一个方便的函数来完成上述步骤。（请注意，LR 测试可以在很少的命令中手动完成！）

---

我对照 R 输出检查了这些结果，接下来给出。

这是R代码：

x <- c(2100,2300,2500,2700,2900,3100,3300,3500,3700,3900,4100,4300)
g <- gl(3, 4)
n <- c(48,42,31,34,31,21,23,23,21,16,17,21)
y <- c(1,2,0,3,8,8,14,17,19,15,17,21)
f <- cbind(y, n-y) ~ x*g
model.matrix(f)  # will be model.frame() for glm()
m1 <- glm(f, family=binomial("probit"))
summary(m1)

以下是完整模型中系数的结果，

Call:
glm(formula = f, family = binomial("probit"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7124  -0.1192   0.1494   0.3036   0.5585  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.892859   2.025096  -1.922   0.0546 .
x            0.000884   0.000829   1.066   0.2863  
g2          -3.212494   2.762155  -1.163   0.2448  
g3          -5.774400   7.553615  -0.764   0.4446  
x:g2         0.001335   0.001017   1.312   0.1894  
x:g3         0.002061   0.002086   0.988   0.3230  

为了比较两个嵌套模型，我使用了以下命令：

m0 <- update(m1, . ~ . -x:g)
anova(m1,m0)

产生以下“偏差表”：

Analysis of Deviance Table

Model 1: cbind(y, n - y) ~ x + g
Model 2: cbind(y, n - y) ~ x * g
  Resid. Df Resid. Dev Df Deviance
1         8     6.4200            
2         6     4.3122  2   2.1078