Uso de Álcool, Cigarro e Maconha - Tabelas Tridimensionais

Não foram utilizados pacotes adicionais na análise. Para a confecção do documento, é necessário o pacote rmarkdown.

Carregamento dos dados:

dados <- array(data = c(911, 3, 44, 2, 538, 43, 456, 279),
               dim = c(2,2,2),
               dimnames = list("Alcohol" = c("Yes","No"),
                               "Cigarette" = c("Yes","No"),
                               "Marijuana" = c("Yes","No")))
ftable(dados, row.vars = c("Marijuana","Alcohol"))

##                   Cigarette Yes  No
## Marijuana Alcohol                  
## Yes       Yes               911  44
##           No                  3   2
## No        Yes               538 456
##           No                 43 279

Verificando as frequências marginais:

addmargins(dados)

## , , Marijuana = Yes
## 
##        Cigarette
## Alcohol Yes No Sum
##     Yes 911 44 955
##     No    3  2   5
##     Sum 914 46 960
## 
## , , Marijuana = No
## 
##        Cigarette
## Alcohol Yes  No  Sum
##     Yes 538 456  994
##     No   43 279  322
##     Sum 581 735 1316
## 
## , , Marijuana = Sum
## 
##        Cigarette
## Alcohol  Yes  No  Sum
##     Yes 1449 500 1949
##     No    46 281  327
##     Sum 1495 781 2276

Verificando as proporções ao longo das linhas:

prop.table(dados, margin = c(1,3))

## , , Marijuana = Yes
## 
##        Cigarette
## Alcohol       Yes        No
##     Yes 0.9539267 0.0460733
##     No  0.6000000 0.4000000
## 
## , , Marijuana = No
## 
##        Cigarette
## Alcohol       Yes        No
##     Yes 0.5412475 0.4587525
##     No  0.1335404 0.8664596

Convertendo para o formato data.frame para usar a função glm e definindo o nível de referência como No para as três variáveis:

dados.df <- as.data.frame(as.table(dados))
dados.df[,-4] <- lapply(dados.df[,-4], relevel, ref = 'No')
dados.df

##   Alcohol Cigarette Marijuana Freq
## 1     Yes       Yes       Yes  911
## 2      No       Yes       Yes    3
## 3     Yes        No       Yes   44
## 4      No        No       Yes    2
## 5     Yes       Yes        No  538
## 6      No       Yes        No   43
## 7     Yes        No        No  456
## 8      No        No        No  279

Modelos Log-Lineares

Serão construídos os seguintes modelos:

Notação	Tipo de associação	Modelo log-linear
(ACM)	Associação conjunta	\(log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ}\)
(AC,AM,CM)	Associação homogênea	\(log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ}\)
(AC,AM), (AC,CM), (AM,CM)	Independência condicional	\(log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ}\)
(A,CM), (C,AM), (M,AC)	Independência conjunta	\(log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY}\)
(A,C,M)	Independência completa	\(log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z\)

mod.ACM <- glm(Freq ~ Cigarette*Marijuana*Alcohol, data = dados.df,
               family = poisson(link = 'log'))
mod.AC.AM.CM <- update(mod.ACM, .~. - Cigarette:Marijuana:Alcohol)
mod.AC.AM <- update(mod.AC.AM.CM, .~. - Cigarette:Marijuana)
mod.AC.CM <- update(mod.AC.AM.CM, .~. - Alcohol:Marijuana)
mod.AM.CM <- update(mod.AC.AM.CM, .~. - Cigarette:Alcohol)
mod.A.CM <- update(mod.AM.CM, .~. - Alcohol:Marijuana)
mod.C.AM <- update(mod.AC.AM, .~. - Cigarette:Alcohol)
mod.M.AC <- update(mod.AC.AM, .~. - Alcohol:Marijuana)
mod.A.C.M <- update(mod.M.AC, .~. - Cigarette:Alcohol)

As contagens estimadas para cada modelo podem ser acessadas usando a função fitted. Assim, pode-se resumir em um único data.frame as contagens estimadas para os modelos ajustados:

estimado <- data.frame(dados.df[,-4],
                       ACM = fitted(mod.ACM),
                       AC.AM.CM = fitted(mod.AC.AM.CM),
                       AC.AM = fitted(mod.AC.AM),
                       AC.CM = fitted(mod.AC.CM),
                       AM.CM = fitted(mod.AM.CM),
                       A.CM = fitted(mod.A.CM),
                       C.AM = fitted(mod.C.AM),
                       M.AC = fitted(mod.M.AC),
                       A.C.M = fitted(mod.A.C.M))
estimado

##   Alcohol Cigarette Marijuana ACM  AC.AM.CM       AC.AM     AC.CM
## 1     Yes       Yes       Yes 911 910.38317 710.0025654 885.87692
## 2      No       Yes       Yes   3   3.61683   0.7033639  28.12308
## 3     Yes        No       Yes  44  44.61683 244.9974346  29.44942
## 4      No        No       Yes   2   1.38317   4.2966361  16.55058
## 5     Yes       Yes        No 538 538.61683 738.9974346 563.12308
## 6      No       Yes        No  43  42.38317  45.2966361  17.87692
## 7     Yes        No        No 456 455.38317 255.0025655 470.55058
## 8      No        No        No 279 279.61683 276.7033639 264.44942
##         AM.CM       A.CM       C.AM      M.AC     A.C.M
## 1 909.2395833 782.682777 627.295694 611.17750 539.98258
## 2   4.7604167 131.317223   3.284271  19.40246  90.59739
## 3  45.7604167  39.391037 327.704306 210.89631 282.09123
## 4   0.2395833   6.608963   1.715729 118.52373  47.32880
## 5 438.8404255 497.525923 652.913005 837.82250 740.22612
## 6 142.1595745  83.474077 211.507030  26.59754 124.19392
## 7 555.1595745 629.400264 341.086995 289.10369 386.70007
## 8 179.8404255 105.599736 110.492970 162.47627  64.87990

Para o modelo de independência condicional entre Álcool e Cigarro, controlado por Maconha (AM,CM), a razão de chances condicional para AC é a mesma nos dois níveis de M:

estimado.AM.CM <- dados
estimado.AM.CM[,,] <- mod.AM.CM$fitted.values
estimado.AM.CM

## , , Marijuana = Yes
## 
##        Cigarette
## Alcohol        Yes         No
##     Yes 909.239583 45.7604167
##     No    4.760417  0.2395833
## 
## , , Marijuana = No
## 
##        Cigarette
## Alcohol      Yes       No
##     Yes 438.8404 555.1596
##     No  142.1596 179.8404

OR <- function(x){x[1,1]*x[2,2]/(x[2,1]*x[1,2])}
apply(estimado.AM.CM,'Marijuana',OR) #AC|M

## Yes  No 
##   1   1

A razão de chances marginal para AC deve desconsiderar M do modelo, ou seja,

marginal <- rowSums(estimado.AM.CM, dims = 2)
OR(marginal)

## [1] 2.749689

Calculando as razões de chances condicionais nos modelos ajustados:

modelos <- list(ACM = mod.ACM,
                AC.AM.CM = mod.AC.AM.CM,
                AC.AM = mod.AC.AM,
                AC.CM = mod.AC.CM,
                AM.CM = mod.AM.CM,
                A.CM = mod.A.CM,
                C.AM = mod.C.AM,
                M.AC = mod.M.AC,
                A.C.M = mod.A.C.M)

ORC <- function(x,condicional){
  est <- dados
  est[,,] <- fitted(x)
  odds.ratio <- apply(est,condicional,OR)
  return(odds.ratio)
}

orc.AC_M <- lapply(modelos,ORC,'Marijuana') #AC|M
orc.AC_M

## $ACM
##       Yes        No 
## 13.803030  7.655141 
## 
## $AC.AM.CM
##      Yes       No 
## 7.803201 7.803201 
## 
## $AC.AM
##    Yes     No 
## 17.703 17.703 
## 
## $AC.CM
##    Yes     No 
## 17.703 17.703 
## 
## $AM.CM
## Yes  No 
##   1   1 
## 
## $A.CM
## Yes  No 
##   1   1 
## 
## $C.AM
## Yes  No 
##   1   1 
## 
## $M.AC
##    Yes     No 
## 17.703 17.703 
## 
## $A.C.M
## Yes  No 
##   1   1

orc.AM_c <- lapply(modelos,ORC,'Cigarette') #AM|C
orc.AM_c

## $ACM
##      Yes       No 
## 24.27076 13.46053 
## 
## $AC.AM.CM
##      Yes       No 
## 19.80658 19.80658 
## 
## $AC.AM
##      Yes       No 
## 61.87324 61.87324 
## 
## $AC.CM
## Yes  No 
##   1   1 
## 
## $AM.CM
##      Yes       No 
## 61.87324 61.87324 
## 
## $A.CM
## Yes  No 
##   1   1 
## 
## $C.AM
##      Yes       No 
## 61.87324 61.87324 
## 
## $M.AC
## Yes  No 
##   1   1 
## 
## $A.C.M
## Yes  No 
##   1   1

orc.CM_A <- lapply(modelos,ORC,'Alcohol') #CM|A
orc.CM_A

## $ACM
##       Yes        No 
## 17.548834  9.732558 
## 
## $AC.AM.CM
##      Yes       No 
## 17.25133 17.25133 
## 
## $AC.AM
## Yes  No 
##   1   1 
## 
## $AC.CM
##     Yes      No 
## 25.1362 25.1362 
## 
## $AM.CM
##     Yes      No 
## 25.1362 25.1362 
## 
## $A.CM
##     Yes      No 
## 25.1362 25.1362 
## 
## $C.AM
## Yes  No 
##   1   1 
## 
## $M.AC
## Yes  No 
##   1   1 
## 
## $A.C.M
## Yes  No 
##   1   1

Calculando as razões de chances marginais nos modelos ajustados:

ORM <- function(x,condicional){
  est <- dados
  est[,,] <- fitted(x)
  est <- aperm(est,condicional)
  odds.ratio<- OR(rowSums(est, dims = 2))
  return(odds.ratio)
}

orm.AC_M <- lapply(modelos,ORM,c(1,2,3)) #AC
orm.AC_M

## $ACM
## [1] 17.703
## 
## $AC.AM.CM
## [1] 17.703
## 
## $AC.AM
## [1] 17.703
## 
## $AC.CM
## [1] 17.703
## 
## $AM.CM
## [1] 2.749689
## 
## $A.CM
## [1] 1
## 
## $C.AM
## [1] 1
## 
## $M.AC
## [1] 17.703
## 
## $A.C.M
## [1] 1

orm.AM_c <- lapply(modelos,ORM,c(1,3,2)) #AM
orm.AM_c

## $ACM
## [1] 61.87324
## 
## $AC.AM.CM
## [1] 61.87324
## 
## $AC.AM
## [1] 61.87324
## 
## $AC.CM
## [1] 5.59619
## 
## $AM.CM
## [1] 61.87324
## 
## $A.CM
## [1] 1
## 
## $C.AM
## [1] 61.87324
## 
## $M.AC
## [1] 1
## 
## $A.C.M
## [1] 1

orm.CM_A <- lapply(modelos,ORM,c(2,3,1)) #CM
orm.CM_A

## $ACM
## [1] 25.1362
## 
## $AC.AM.CM
## [1] 25.1362
## 
## $AC.AM
## [1] 1.932727
## 
## $AC.CM
## [1] 25.1362
## 
## $AM.CM
## [1] 25.1362
## 
## $A.CM
## [1] 25.1362
## 
## $C.AM
## [1] 1
## 
## $M.AC
## [1] 1
## 
## $A.C.M
## [1] 1

Os valores estimados dos parâmetros podem ser usados para obter as razões de chances condicionais. Por exemplo, para o modelo (AM,CM), tem-se:

summary(mod.AM.CM)

## 
## Call:
## glm(formula = Freq ~ Cigarette + Marijuana + Alcohol + Cigarette:Marijuana + 
##     Marijuana:Alcohol, family = poisson(link = "log"), data = dados.df)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  0.0584  -0.8663  -0.2619   2.2287   4.5702  -9.7716  -4.3441   6.8353  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                5.19207    0.06088  85.285  < 2e-16 ***
## CigaretteYes              -0.23512    0.05551  -4.235 2.28e-05 ***
## MarijuanaYes              -6.62092    0.47370 -13.977  < 2e-16 ***
## AlcoholYes                 1.12719    0.06412  17.579  < 2e-16 ***
## CigaretteYes:MarijuanaYes  3.22431    0.16098  20.029  < 2e-16 ***
## MarijuanaYes:AlcoholYes    4.12509    0.45294   9.107  < 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: 2851.46  on 7  degrees of freedom
## Residual deviance:  187.75  on 2  degrees of freedom
## AIC: 248.8
## 
## Number of Fisher Scoring iterations: 5

coeficientes <- mod.AM.CM$coefficients
coeficientes

##               (Intercept)              CigaretteYes 
##                 5.1920699                -0.2351197 
##              MarijuanaYes                AlcoholYes 
##                -6.6209239                 1.1271857 
## CigaretteYes:MarijuanaYes   MarijuanaYes:AlcoholYes 
##                 3.2243089                 4.1250878

A razão de chances condicional para AC no modelo (AM,CM) é \(e^{0} = 1\).

A razão de chances condicional para AM no modelo (AM,CM) é \(e^{4.125} = 61.873\).

A razão de chances condicional para CM no modelo (AM,CM) é \(e^{3.224} = 25.136\).

Qualidade do ajuste

As estatísticas \(\chi^2\) e \(G^2\) podem ser usadas para comparar os valores observados com os ajustados (para grandes amostras).

Uma razão elevada para \(\chi^2/gl\) e \(G^2/gl\), com \(gl\): número de graus de liberdade (diferença entre o número de células e o número de parâmetros no modelo), indica rejeição de \(H_0\): Modelo log-linear é adequado.

Para os modelos ajustados, tem-se:

GoF <- function(x){
  G2 <- round(x$deviance, 1)
  X2 <- round(sum(resid(x, type = "pearson")^2), 1)
  df <- round(x$df.residual, 1)
  p.valor <- round(1 - pchisq(G2,df), 1)
  AIC <- round(x$aic, 1)
  return(c(G2 = G2, X2 = X2, gl = df, p.valor.G2 = p.valor, AIC = AIC))
}

QdA <- lapply(modelos, GoF)
t(rbind.data.frame(QdA))

##              G2     X2 gl p.valor.G2    AIC
## ACM         0.0    0.0  0        1.0   65.0
## AC.AM.CM    0.4    0.4  1        0.5   63.4
## AC.AM     497.4  443.8  2        0.0  558.4
## AC.CM      92.0   80.8  2        0.0  153.1
## AM.CM     187.8  177.6  2        0.0  248.8
## A.CM      534.2  505.6  3        0.0  593.3
## C.AM      939.6  824.2  3        0.0  998.6
## M.AC      843.8  704.9  3        0.0  902.9
## A.C.M    1286.0 1411.4  4        0.0 1343.1

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/.

Uso de Álcool, Cigarro e Maconha - Tabelas Tridimensionais

Modelos Log-Lineares - Análise de Dados Categorizados - Prova 3

Allan Vieira - 14/0128492

Caio Balena - 13/0104566

Frederico de Lucca - 15/0072139

José Cezário - 15/0072414

Junho de 2018

Modelos Log-Lineares

Qualidade do ajuste

Referências Bibliográficas