Codificación de una interacción entre un predictor nominal y continuo de regresión logística en MATLAB

Entonces nuestros datos están estructurados de la siguiente manera:

Tenemos participantes, cada participante se puede clasificar en 3 grupos ( G ), y para cada participante tenemos muestras de una variable continua. Y estamos tratando de predecir valores que son 0 o 1.$M$ $\in {A,B,C}$$N$

¿Cómo usaríamos matlab para probar una interacción entre la variable continua y la variable categórica al predecir estos valores?

La forma más fácil, IMO, es construir la matriz de diseño usted mismo, ya que glmfitacepta una matriz de valores brutos (observados) o una matriz de diseño. Codificar un término de interacción no es tan difícil una vez que escribió el modelo completo. Digamos que tenemos dos predictores, (continuo) (categórico, con tres niveles desordenados, digamos ). Usando la notación de Wilkinson, escribiríamos este modelo como , descuidando el lado izquierdo (para un resultado binomial, usaríamos una función de enlace logit). Solo necesitamos dos vectores ficticios para codificar los niveles (como presente / ausente para una observación particular), por lo que tendremos 5 coeficientes de regresión, más un término de intercepción. Esto se puede resumir como$x$$g$$g={1,2,3}$y ~ x + g + x:gg

donde representa una matriz indicadora que codifica el nivel de .$\mathbb{I}$$g$

En Matlab, usando el ejemplo en línea, haría lo siguiente:

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');

No incluí una columna de unos para la intercepción ya que está incluida por defecto. La matriz de diseño se parece a

    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

y puede ver que los términos de interacción están codificados como el producto de xcon la columna correspondiente de g(g = 2 y g = 3, ya que no necesitamos el primer nivel).

Los resultados se dan a continuación, como coeficientes, errores estándar, estadística y valor p (de la statsestructura):

   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

Ahora, se puede probar la interacción calculando la diferencia en la desviación del modelo completo anterior y un modelo reducido (omitiendo el término de interacción, es decir, las dos últimas columnas de la matriz de diseño). Esto se puede hacer manualmente o utilizando la lratiotestfunción que proporciona la prueba de hipótesis de la razón de verosimilitud. La desviación para el modelo completo es 4.3122 ( dev), mientras que para el modelo sin interacción es 6.4200 (utilicé glmfit(X(:,1:3), [y n], 'binomial', 'link', 'probit');), y la prueba LR asociada tiene dos grados de libertad (la diferencia en el número de parámetros entre los dos modelos). Como la desviación escalada es solo dos veces la probabilidad logarítmica de GLM, podemos usar

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

donde la estadística se distribuye como a con 2 df (el valor crítico es entonces 5.9915, ver ). El resultado indica un resultado no significativo: no podemos concluir la existencia de una interacción entre y en la muestra observada.$\chi^2$chi2inv(0.95, 2)xg

Supongo que puede concluir los pasos anteriores en una función conveniente de su elección. (¡Tenga en cuenta que la prueba LR puede hacerse a mano en muy pocos comandos!)

Verifiqué esos resultados con la salida R, que se proporciona a continuación.

Aquí está el código 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)

Aquí están los resultados, para los coeficientes en el modelo completo,

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  

Para la comparación de los dos modelos anidados, utilicé los siguientes comandos:

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

que produce la siguiente "tabla de desviaciones":

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