在 MATLAB 中為邏輯回歸編碼名義變量和連續預測變量之間的交互
所以我們的數據結構如下：
我們有參與者，每個參與者可分為 3 組（G )，對於每個參與者，我們有連續變量的樣本。我們正在嘗試預測 0 或 1 的值。
在預測這些值時，我們如何使用 matlab 來測試連續變量和分類變量之間的交互作用？
IMO 最簡單的方法是自己構建設計矩陣，因為它
glmfit接受原始（觀察到的）值矩陣或設計矩陣。一旦你編寫了完整的模型，編寫一個交互項就不是那麼困難了。假設我們有兩個預測變量，（連續）和（分類，具有三個無序級別，例如）。使用 Wilkinson 的符號，我們將把這個模型寫成y ~ x + g + x:g，忽略左側（對於二項式結果，我們將使用 logit 鏈接函數）。我們只需要兩個虛擬向量來編碼g級別（特定觀察的存在/不存在），所以我們將有 5 個回歸係數，加上一個截距項。這可以概括為在哪裡代表一個指標矩陣，編碼的水平.
在 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)其中統計量分佈為使用 2 df（臨界值為 5.9915，請參閱
chi2inv(0.95, 2)）。輸出表明一個不顯著的結果：我們不能斷定觀察到的樣本x之間存在相互作用。g我想你可以用你選擇的一個方便的函數來完成上述步驟。（請注意，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