Carregamento de pacotes:

Neste exemplo é utilizado o pacote: vcdExtra. Para a confecção do documento, é necessário o pacote rmarkdown.

# vcdExtra -- para o Teste de Mantel Haenzel
if (!require('vcdExtra')) install.packages('vcdExtra'); library('vcdExtra')

Carregamento dos dados:

dados <- array(data = c(1339, 260, 88, 300, 55, 22),
               dim = c(3, 2),
               dimnames = list('Race' = c('White', 'Black', 'Others'),
                               'Belief' = c('Yes', 'No')))
dados
##         Belief
## Race      Yes  No
##   White  1339 300
##   Black   260  55
##   Others   88  22

Verificando as frequências marginais:

addmargins(dados)
##         Belief
## Race      Yes  No  Sum
##   White  1339 300 1639
##   Black   260  55  315
##   Others   88  22  110
##   Sum    1687 377 2064

Verificando as proporções ao longo das linhas:

prop.table(dados, margin = 1)
##         Belief
## Race           Yes        No
##   White  0.8169616 0.1830384
##   Black  0.8253968 0.1746032
##   Others 0.8000000 0.2000000

Algumas estatísticas para a tabela Bidimensional:

vcdExtra::CMHtest(dados)
## Cochran-Mantel-Haenszel Statistics for Race by Belief 
## 
##                  AltHypothesis    Chisq Df    Prob
## cor        Nonzero correlation 0.017785  1 0.89391
## rmeans  Row mean scores differ 0.359901  2 0.83531
## cmeans  Col mean scores differ 0.017785  1 0.89391
## general    General association 0.359901  2 0.83531

Convertendo para o formato data.frame para usar a função glm e definindo o último nível como nível de referência:

dados.df <- as.data.frame(as.table(dados))
dados.df$Race <- relevel(dados.df$Race, ref = 'Others')
dados.df$Belief <- relevel(dados.df$Belief, ref = 'No')
dados.df
##     Race Belief Freq
## 1  White    Yes 1339
## 2  Black    Yes  260
## 3 Others    Yes   88
## 4  White     No  300
## 5  Black     No   55
## 6 Others     No   22

Modelo de Independência

Construindo o modelo log-linear de independência da forma \(log(\mu_{ij}) = \lambda + \lambda_i^X + \lambda_j^Y \quad i = 1,2,3 \quad j = 1,2\), com \(\lambda\): efeito geral, \(\lambda_i^X\): efeito marginal da raça (1: White, 2: Black, 3: Others) e \(\lambda_j^Y\): efeito marginal da crença (1: Yes, 2: No), tem-se:

mod.ind <- glm(Freq ~ Race + Belief, data = dados.df, family = poisson(link = 'log'))
summary(mod.ind)
## 
## Call:
## glm(formula = Freq ~ Race + Belief, family = poisson(link = "log"), 
##     data = dados.df)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6  
## -0.01717   0.15781  -0.20194   0.03631  -0.33688   0.41917  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.00032    0.10611   28.28   <2e-16 ***
## RaceWhite    2.70136    0.09849   27.43   <2e-16 ***
## RaceBlack    1.05209    0.11075    9.50   <2e-16 ***
## BeliefYes    1.49846    0.05697   26.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2849.21758  on 5  degrees of freedom
## Residual deviance:    0.35649  on 2  degrees of freedom
## AIC: 49.437
## 
## Number of Fisher Scoring iterations: 3
coeficientes <- summary(mod.ind)$coef

O R utiliza restrições do tipo set-to-zero (codificação dummy), definindo por construção o nível de referência igual a zero e portanto, não mostrando seu valor junto aos coeficientes estimados.

Para o modelo log-linear de independência, com restrição set-to-zero no último nível (referência), os coeficientes estimados são \(\lambda^X_1\) = 2.7013612, \(\lambda^X_2\) = 1.0520923, \(\lambda^X_3\) = 0 (referência), \(\lambda^Y_1\) = 1.4984619 e \(\lambda^Y_2\) = 0 (referência).

Assim, a chance estimada de crença em vida após a morte é igual a \(e^{\lambda^Y_1} = e^{1.498} = 4.475\).

Como foi usado um modelo de independência, a razão de chance é igual a 1 para toda sub-tabela 2x2 das contagens estimadas:

estimado <- dados
estimado[,] <- mod.ind$fitted.values
estimado
##         Belief
## Race            Yes        No
##   White  1339.62839 299.37161
##   Black   257.46366  57.53634
##   Others   89.90795  20.09205
estimado[1,1]*estimado[2,2]/(estimado[2,1]*estimado[1,2])
## [1] 1
estimado[1,1]*estimado[3,2]/(estimado[3,1]*estimado[1,2])
## [1] 1
estimado[2,1]*estimado[3,2]/(estimado[3,1]*estimado[2,2])
## [1] 1

Qualidade do ajuste:

A estatística \(G^2\) é mostrada diretamente no resultado do modelo ajustado:

G2 <- mod.ind$deviance # estatística G2
G2
## [1] 0.3564851
G2.df <- mod.ind$df.residual # graus de liberdade
G2.df
## [1] 2
G2.pv <- 1 - pchisq(G2, G2.df) # p-valor
G2.pv
## [1] 0.8367395
G2/G2.df
## [1] 0.1782425

A estatística \(\chi^2\) pode ser obtida fazendo:

X2 <- sum(resid(mod.ind, type = "pearson")^2) # estatística X2
X2
## [1] 0.3600752
X2.df <- mod.ind$df.residual # graus de liberdade
X2.df
## [1] 2
X2.pv <- 1 - pchisq(X2, X2.df) # p-valor
X2.pv
## [1] 0.8352388
X2/X2.df
## [1] 0.1800376

Modelo Saturado

Construindo o modelo log-linear saturado da forma \(log(\mu_{ij}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_{ij}^{XY} \quad i = 1,2,3 \quad j = 1,2\), com \(\lambda\): efeito geral, \(\lambda_i^X\): efeito marginal da raça (1: White, 2: Black, 3: Others), \(\lambda_j^Y\): efeito marginal da crença (1: Yes, 2: No) e \(\lambda_{ij}^{XY}\): parâmetros de associação que refletem o desvio de independência, tem-se:

mod.sat <- glm(Freq ~ Race*Belief, data = dados.df, family = poisson(link = 'log'))
summary(mod.sat)
## 
## Call:
## glm(formula = Freq ~ Race * Belief, family = poisson(link = "log"), 
##     data = dados.df)
## 
## Deviance Residuals: 
## [1]  0  0  0  0  0  0
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           3.0910     0.2132  14.498  < 2e-16 ***
## RaceWhite             2.6127     0.2209  11.829  < 2e-16 ***
## RaceBlack             0.9163     0.2523   3.632 0.000281 ***
## BeliefYes             1.3863     0.2384   5.816 6.03e-09 ***
## RaceWhite:BeliefYes   0.1096     0.2468   0.444 0.656946    
## RaceBlack:BeliefYes   0.1671     0.2808   0.595 0.551889    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  2.8492e+03  on 5  degrees of freedom
## Residual deviance: -9.4813e-14  on 0  degrees of freedom
## AIC: 53.081
## 
## Number of Fisher Scoring iterations: 3
coeficientes <- summary(mod.sat)$coef

A razão de chances estimada entre crença e raça é igual a \(e^{\lambda^{XY}_{11}} = e^{0.1096} = 1.1158\) para brancos e outros e \(e^{\lambda^{XY}_{21}} = e^{0.1671} = 1.1818\).

Assim, a razão de chances estimada entre crença e raça é igual a \(e^{\lambda^{XY}_{11} - \lambda^{XY}_{21}} = e^{-0.0575} = 0.9442\) para brancos e negros.

Para o modelo saturado, as estatísticas \(\chi^2\) e \(G^2\) são nulas, com número de graus de liberdade igual a zero.

Referências Bibliográficas

AGRESTI, ALAN. An Introduction to Categorical Data Analysis. John Wiley & Sons, second edition, 2007.

NOTAS DE AULA. Análise de Dados Categorizados. Curso de Graduação em Estatística, UnB, 2018.

R CORE TEAM. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.