Codificando uma interação entre um preditor nominal e um preditor contínuo para regressão logística no MATLAB

Portanto, nossos dados estão estruturados da seguinte forma:

Temos participantes, cada participante pode ser categorizado em 3 grupos ( G \ em {A, B, C} ) e, para cada participante, temos N amostras de uma variável contínua. E estamos tentando prever valores que são 0 ou 1.$M$ $\in {A,B,C}$$N$

Como usaríamos o matlab para testar uma interação entre a variável contínua e a variável categórica na previsão desses valores?

A maneira mais fácil, IMO, é construir você mesmo a matriz de design, pois glmfitaceita uma matriz de valores brutos (observados) ou uma matriz de design. Codificar um termo de interação não é tão difícil assim que você escreve o modelo completo. Digamos que temos dois preditores, (contínuo) (categórico, com três níveis não ordenados, digamos ). Usando a notação de Wilkinson, escreveríamos esse modelo como negligenciando o lado esquerdo (para um resultado binomial, usaríamos uma função de logit link). Precisamos apenas de dois vetores fictícios para codificar os níveis (como presente / ausente para uma observação em particular); portanto, teremos 5 coeficientes de regressão, além de um termo de interceptação. Isso pode ser resumido como$x$$g$$g={1,2,3}$y ~ x + g + x:gg

onde representa uma matriz indicadora que codifica o nível de .$\mathbb{I}$$g$

No Matlab, usando o exemplo online, eu faria o seguinte:

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

Não incluí uma coluna de uns para a interceptação, pois ela é incluída por padrão. A matriz de design parece

    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

e você pode ver que os termos de interação são apenas codificados como produto da xcoluna correspondente de g(g = 2 eg = 3, pois não precisamos do primeiro nível).

Os resultados são apresentados abaixo, como coeficientes, erros padrão, estatística e valor-p (da statsestrutura):

   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

Agora, o teste da interação pode ser feito calculando a diferença de desvio do modelo completo acima e um modelo reduzido (omitindo o termo de interação, que são as duas últimas colunas da matriz de design). Isso pode ser feito manualmente ou usando a lratiotestfunção que fornece o teste de hipótese da razão de verossimilhança. O desvio para o modelo completo é 4.3122 ( dev), enquanto para o modelo sem interação é 6.4200 (eu usei glmfit(X(:,1:3), [y n], 'binomial', 'link', 'probit');), e o teste LR associado tem dois graus de liberdade (a diferença no número de parâmetros entre os dois modelos). Como o desvio em escala é apenas duas vezes a probabilidade de log para GLMs, podemos usar

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

onde a estatística é distribuída como a com 2 df (o valor crítico é então 5,9915, consulte ). O resultado indica um resultado não significativo: Não podemos concluir a existência de uma interação entre e na amostra observada.$\chi^2$chi2inv(0.95, 2)xg

Eu acho que você pode concluir as etapas acima em uma função conveniente de sua escolha. (Observe que o teste LR pode ser feito manualmente em muito poucos comandos!)

Eu verifiquei esses resultados com relação à saída R, que é dada a seguir.

Aqui está o 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)

Aqui estão os resultados, para os coeficientes no 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 a comparação dos dois modelos aninhados, usei os seguintes comandos:

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

que produz a seguinte "tabela de desvio":

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