Modelo Linear Generalizado - Distribuição Multinomial

Valeska Andreozzi (CEAUL) 1

Modelo Linear Generalizado

Distribuicao Multinomial

Problema 1

Os dados a seguir corresponde a 3417 mulheres em idade fertil categorizadas por faixa etariae metodo contraceptivo.

> cuse <- read.table("brazil.dat", header = F)

> names(cuse) <- c("age", "met", "n")

> cuse

age met n

1 15-19 ster 2

2 15-19 other 75

3 15-19 no 90

4 20-24 ster 32

5 20-24 other 223

6 20-24 no 280

7 25-29 ster 147

8 25-29 other 296

9 25-29 no 237

10 30-34 ster 298

11 30-34 other 197

12 30-34 no 207

13 35-39 ster 276

14 35-39 other 96

15 35-39 no 196

16 40-44 ster 202

17 40-44 other 71

18 40-44 no 174


> summary(cuse)

age met n

15-19:3 no :6 Min. : 2.0

20-24:3 other:6 1st Qu.: 91.5

25-29:3 ster :6 Median :196.5

30-34:3 Mean :172.2

35-39:3 3rd Qu.:233.5

40-44:3 Max. :298.0

A tabela a seguir mostra o quao associadas estao as duas variaveis (idade e metodo decontracepcao)

> library(gregmisc)

> indiv <- data.frame(rep(cuse$age, cuse$n), rep(cuse$met,

+ cuse$n))

> tab <- CrossTable(indiv[, 1], indiv[, 2], prop.r = F,

+ prop.c = F, prop.t = F, prop.chisq = F, chisq = T,

+ dnn = c("age", "metodo"))


Cell Contents

|-------------------------|

| N |

| N / Row Total |

|-------------------------|

Total Observations in Table: 3099

| metodo

age | no | other | ster | Row Total |

-------------|-----------|-----------|-----------|-----------|

15-19 | 90 | 75 | 2 | 167 |

| 0.539 | 0.449 | 0.012 | 0.054 |

-------------|-----------|-----------|-----------|-----------|

20-24 | 280 | 223 | 32 | 535 |

| 0.523 | 0.417 | 0.060 | 0.173 |

-------------|-----------|-----------|-----------|-----------|

25-29 | 237 | 296 | 147 | 680 |

| 0.349 | 0.435 | 0.216 | 0.219 |

-------------|-----------|-----------|-----------|-----------|

30-34 | 207 | 197 | 298 | 702 |

| 0.295 | 0.281 | 0.425 | 0.227 |

-------------|-----------|-----------|-----------|-----------|

35-39 | 196 | 96 | 276 | 568 |

| 0.345 | 0.169 | 0.486 | 0.183 |

-------------|-----------|-----------|-----------|-----------|

40-44 | 174 | 71 | 202 | 447 |

| 0.389 | 0.159 | 0.452 | 0.144 |

-------------|-----------|-----------|-----------|-----------|

Column Total | 1184 | 958 | 957 | 3099 |

-------------|-----------|-----------|-----------|-----------|

Statistics for All Table Factors

Pearson’s Chi-squared test

------------------------------------------------------------

Chi^2 = 484.7404 d.f. = 10 p = 8.032588e-98


Atraves do teste χ2 testamos a probabilidade conjunta da variavel Y (metodo contracep-tivo) e X (idade), isto e, testamos a independencia.

Mas como podemos investigar se idade (X) e preditor da variavel resposta metodo con-traceptivo (Y )? Agora estamos interessados na distribuicao condicional da variavel respostaY dado o preditor X.

Modelos Multinomiais

Distribuicao multinomial

Considere Yi uma variavel resposta que pode assumir valores discretos, nomeadamente,1, 2, · · · , J . Assumindo que a nossa variavel no exemplo anterior e metodo contraceptivotemos:

pij = Pr(Yi = j)

como sendo a probabilidade de uma mulher pertencer a algumas das categorias de utilizacaode contraceptivo (esterilizacao, outro metodo, nenhum metodo). Caracterısticas:

• as categorias sao mutuamente exclusivas

• para cada indivıduo a soma das probabilidades e igual a 1 (∑j

j=1pij = 1)

Para dados agrupados, como o nosso exemplo, temos uma variavel ni que denota onumero de casos do grupo i e seja Yij o numero de casos do grupo i na categoria j. O ındicei representa o grupo idade e j as categorias de metodo contraceptivo. Logo, tem-se que,y11, y12, y13, representam o numero de mulheres entre 15 e 19 anos que sao esterilizadas, queutilizam outro metodo de contracepcao ou que nao usam nenhum metodo, respectivamente.

Para dados individuais temos ni = 1 e y11 = 1 ou y11 = 0 se uma mulher e esterilizadaou nao, respectivamente.

A distribuicao de probabilidade para a variavel Yij dado o total de casos ni e dado por:

Pr(Yi1 = yi1, · · · , Yij = yij) =

(ni

yi1, · · · , yij

)pyi1i1 × · · · × p

yijij (1)

Quando j = 2 temos a distribuicao binomial.


Modelo Logit Multinomial

Como podemos investigar se idade (X) e preditor da variavel resposta metodo contraceptivo(Y )?

Podemos assumir que uma das categorias da variavel resposta e referencia e calcular achance (odds) das outras categorias em relacao a esta categoria referencia. No nosso exemplo,vamos assumir a categoria nao utiliza contraceptivo” como sendo a categoria referencia.Desta forma terıamos a chance de “esterilizacao” comparado com a chance de “nao utilizacaode metodo contraceptivo” e a chance de “outro metodo” comparado com a chance de “naoutilizacao de metodo contraceptivo”.

Veja o grafico do logit (logarıtmo da chance) empırico dos dados:

> oddsster <- log(tab$prop.row[, 3]/tab$prop.row[,

+ 1])

> oddsster

15-19 20-24 25-29 30-34 35-39

-3.8066625 -2.1690537 -0.4776276 0.3643747 0.3422862

40-44

0.1492124

> oddsother <- log(tab$prop.row[, 2]/tab$prop.row[,

+ 1])

> oddsother

15-19 20-24 25-29 30-34 35-39

-0.18232156 -0.22761783 0.22229931 -0.04951506 -0.71376647

40-44

-0.89637542

> plot(1:6, oddsster, ylim = range(oddsster, oddsother),

+ pch = 19, col = "blue", axes = F, xlab = "idade",

+ ylab = "logit")

> points(1:6, oddsother, pch = 19, col = "red")

> lines(lowess(1:6, oddsster, f = 0.3), col = "blue")

> lines(lowess(1:6, oddsother, f = 0.3), col = "red")


> box()

> axis(2)

> axis(1, at = 1:6, labels = levels(cuse$age))

> legend(4, -3, legend = c("esterilizac~ao", "outro metodo"),

+ bty = "n", pch = 19, col = c("blue", "red"))

idade

logi

t

−3

−2

−1

0

15−19 20−24 25−29 30−34 35−39 40−44

esterilizaçãooutro método

O grafico do logit empırico sugere que o logaritmo da chance das mulheres esterilizadasaumenta rapido com a idade, atinge o maximo aos 30-34 anos e declina lentamente a partirdessa faixa etaria. O logarıtmo da chance de outro metodo atinge o maximo aos 25-29 edecresce rapidamente a seguir.

Seja 1 = esterilizacao, 2 = outro metodo, 3 = nenhum metodo. O modelo logit multino-mial e igual a

log

(pi1pi3

)= β0 + β1agei

log

(pi2pi3

)= α0 + α1agei (2)


Para a variavel resposta multinomial temos j − 1 regressoes logısticas. Genericamente,podemos definir o modelo logıstico multinomial da seguinte forma:

log

(pijpiJ

)= αj + x′

iβj (3)

em que j = 1, . . . , J − 1, xi e o vetor de covariaveis, αj sao os interceptos e βj e o vetor decoeficientes de regressao.

Na formula (2) estamos comparando categoria 1:3 e 2:3. Se quisermos comparar categorias1:2 basta fazer

log

(pi1pi2

)= log

(pi1pi3

)− log

(pi2pi3

)

O grafico dos logits sugere um relacao quadratica com a idade. Para contemplar estacaracterıstica observada nos dados vamos criar uma variavel idade que e igual ao pontomedio de cada faixa etaria.

> cuse$idade <- rep(seq(17, 42, length = 6), each = 3)

E vamos ajustar seguinte modelo:

log

(pi1pi3

)= β0 + β1agei + β2age

2

i

log

(pi2pi3

)= α0 + α1agei + α2age

2

i

> library(nnet)

> mod1 <- multinom(met ~ idade + I(idade^2), data = cuse,

+ weights = n)

# weights: 12 (6 variable)

initial value 3404.599483

iter 10 value 3112.762332

final value 3112.745366

converged

> mod1


Call:

multinom(formula = met ~ idade + I(idade^2), data = cuse, weights = n)

Coefficients:

(Intercept) idade I(idade^2)

other -3.042142 0.2340215 -0.004432943

ster -16.071596 0.9166541 -0.012667215

Residual Deviance: 6225.491

AIC: 6237.491

> round(unique(mod1$fit), 3)

no other ster

1 0.580 0.411 0.009

4 0.478 0.460 0.062

7 0.382 0.399 0.219

10 0.308 0.281 0.412

13 0.304 0.193 0.503

16 0.409 0.146 0.445

Podemos comparar o mod1 com um modelo mais simples, retirando o termo ao quadradoda idade

> mod2 <- multinom(met ~ idade, data = cuse, weights = n)

# weights: 9 (4 variable)

initial value 3404.599483

final value 3184.024055

converged

> mod2

Call:

multinom(formula = met ~ idade, data = cuse, weights = n)


Coefficients:

(Intercept) idade

other 0.7614718 -0.03357687

ster -3.4062103 0.09929226

Residual Deviance: 6368.048

AIC: 6376.048

> anova(mod1, mod2)

Likelihood ratio tests of Multinomial Models

Response: met

Model Resid. df Resid. Dev Test Df

1 idade 32 6368.048

2 idade + I(idade^2) 30 6225.491 1 vs 2 2

LR stat. Pr(Chi)

1

2 142.5574 0

E se quisermos as probabilidades? Como calcular?

log

(pi1pi3

)= β0 + β1agei + β2age

2

i = η1

log

(pi2pi3

)= α0 + α1agei + α2age

2

i = η2

Consequentemente temos:

pi1pi3

= exp(β0 + β1agei + β2age2

i ) = exp(η1) (4)

pi2pi3

= exp(α0 + α1agei + α2age2

i ) = exp(η2) (5)

E alem disso temos que:pi1 + pi2 + pi3 = 1 (6)


Substituindo (4) e (5) em (6) temos que:

pi1 =exp(η1)

1 + exp(η1) + exp(η2)

pi2 =exp(η2)

1 + exp(η1) + exp(η2)

pi3 =1

1 + exp(η1) + exp(η2)

A figura a seguir mostra o grafico das probabilidades estimado pelo modelo mod1

> predicao <- predict(mod1, cuse[, c(2, 4)], type = "probs")

> prob <- round(as.vector(t(unique(predicao)[, c(3,

+ 2, 1)])), 2)

> pred <- cbind(cuse[, c(1:2, 4)], prob)


> pred

age met idade prob

1 15-19 ster 17 0.01

2 15-19 other 17 0.41

3 15-19 no 17 0.58

4 20-24 ster 22 0.06

5 20-24 other 22 0.46

6 20-24 no 22 0.48

7 25-29 ster 27 0.22

8 25-29 other 27 0.40

9 25-29 no 27 0.38

10 30-34 ster 32 0.41

11 30-34 other 32 0.28

12 30-34 no 32 0.31

13 35-39 ster 37 0.50

14 35-39 other 37 0.19

15 35-39 no 37 0.30

16 40-44 ster 42 0.44

17 40-44 other 42 0.15

18 40-44 no 42 0.41

> matplot(y = predicao, x = pred$idade, type = "l",

+ col = c(1, 2, 4), lty = 1, lwd = 1.5, ylab = "Probabilidade",

+ xlab = "", axes = F)

> axis(2)

> box()

> axis(1, at = unique(cuse$idade), labels = levels(cuse$age))

> legend(30, 0.15, c("no", "other", "ster"), col = c(1,

+ 2, 4), lty = 1, lwd = 1.5)


Pro

babi

lidad

e

0.0

0.1

0.2

0.3

0.4

0.5

0.6

15−19 20−24 25−29 30−34 35−39 40−44

nootherster

Estimacao

O modelo multinomial acima descrito e estimado atraves do metodo da maxima verossimi-lhanca. O metodo consiste em maximizar a verossimilhanca em (2) com probabilidades pijvisto como funcao dos parametros de regressao β′s.

Modelo para variavel ordinal

Em muitas situacoes da vida pratica temos uma variavel resposta ordinal ao inves de nominal(como no exemplo anterior). A ordem presente na variavel resposta Y pode acontecer porela representar uma escala, como por exemplo, severidade de dor (nenhuma, leve, moderada,grave). Em outros casos a variavel resposta ordinal e construıda atraves da especificacaode hierarquias (exemplo: anos de escolaridade). Existem varios modelos para modelar umavariavel resposta ordinal. Neste curso veremos dois tipos: modelo de chances proporcionais(proportional odds) tambem conhecido como modelo cumulativo (cumulative logit model) emodelo de razao contınua (continuation ratio model).


Para efeito de ilustracao, vamos considerar o seguinte exemplo: investigar fatores asso-ciados a satisfacao das condicoes da casa para residentes de Copenhagen. Os dados estaodisponıveis no formato agrupado.

> library(rms)

> copen <- read.table("copen.dat")

> head(copen)

housing influence contact satisfaction n

1 tower low low low 21

2 tower low low medium 21

3 tower low low high 28

4 tower low high low 14

5 tower low high medium 19

6 tower low high high 37

As variaveis sao as seguintes:

• housing tipo de casa (tower blocks, apartments, atrium houses and terraced houses),

• influence sentimento de influencia no gerenciamento do apartamento (low, medium,high),

• contact grau de contato com os vizinhos (low, high), and

• satisfaction satisfacao com as condicoes da casa (low, medium, high)

• n numero de residentes nas diversas combinacoes das variaveis acima descritas

Modelo de chances proporcionais

Para uma variavel resposta categorica ordinal Y com 1, 2, . . . , k nıveis temos

πij = Pr(Yi = j)

a probabilidade da variavel resposta pertencer a categoria j = 1, 2, . . . , k e


γij = Pr(Yi ≤ j) = πi1 + πi2 + · · ·+ πij

a probabilidade cumulativa correspondente de Yi. O modelo de chances proporcionais paraYi e descrito por:

Pr(Yi ≤ j|Xi) =1

1 + exp[−(αj + x′

iβ)](7)

Podemos reescrever o modelo (7) da seguinte forma:

logit(γij) = lnγij

1− γij= αj + x′

iβ (8)

O modelo de chances proporcionais compara a probabilidade de uma resposta Y ser igualou menor que uma determinada categoria com a probabilidade da resposta ser maior queesta categoria. Alem disso, o modelo e composto por k − 1 equacoes lineares paralelas.Para um nıvel j fixo ou para uma categoria fixa, o modelo torna-se uma regressao logısticasimples para o evento Y ≤ j. Nota-se que o modelo (7) contem k − 1 + p parametros (pe o numero de covariaveis no modelo e k e o numero de categorias de Y ). Cada equacaopossui um intercepto que atende a seguinte condicao α1 < α2 < . . . < αk. Os parametrosβ nao depende da categoria j. Esse modelo fornece uma unica estimativa para a razao dechances para todas as categorias comparadas, que pode ser obtida atraves do calculo doexponencial do coeficiente β. Essa premissa e bastante adequada em termos da facilidade deinterpretacao e de parcimonia do modelo.

Vejamos o modelo de chances porporcionais estimado para os dados de satisfacao. Paratal temos que transformar os dados agrupados em dados individuais:

> coindiv <- data.frame(housing = rep(copen$housing,

+ copen$n), influence = rep(copen$influence,

+ copen$n), contact = rep(copen$contact, copen$n),

+ satisfaction = rep(copen$satisfaction, copen$n))

A seguir criamos uma variavel resposta com labels numericos para garantir a hierarquiada variavel resposta ordinal

> coindiv$satis2 <- 1

> coindiv$satis2[which(coindiv$satisfaction == "medium")] <- 2

> coindiv$satis2[which(coindiv$satisfaction == "high")] <- 3

> coindiv$satis2 <- factor(coindiv$satis2)

> library(VGAM)

> fit1 <- vglm(satis2 ~ housing, cumulative(parallel = TRUE),

+ coindiv)

> summary(fit1)


Call:

vglm(formula = satis2 ~ housing, family = cumulative(parallel = TRUE),

data = coindiv)

Pearson Residuals:

Min 1Q Median 3Q Max

logit(P[Y<=1]) -1.3612 -0.9960 -0.35744 1.3221 1.6806

logit(P[Y<=2]) -1.6285 -1.0954 0.40890 1.0098 1.3922

Coefficients:

Value Std. Error t value

(Intercept):1 -0.64981 0.071302 -9.1134

(Intercept):2 0.47151 0.070384 6.6991

housingatrium -0.18811 0.137169 -1.3713

housingterraced 0.58502 0.131241 4.4576

housingtower -0.46842 0.115644 -4.0506

Number of linear predictors: 2

Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2])

Dispersion Parameter for cumulative family: 1

Residual Deviance: 3594.803 on 3357 degrees of freedom

Log-likelihood: -1797.401 on 3357 degrees of freedom

Number of Iterations: 3

> fit2 <- vglm(satis2 ~ housing + influence, cumulative(parallel = TRUE),

+ coindiv)

> summary(fit2)

Call:

vglm(formula = satis2 ~ housing + influence, family = cumulative(parallel = TRUE),

data = coindiv)

Pearson Residuals:



logit(P[Y<=1]) -1.5618 -0.87983 -0.30840 1.02639 2.5220

logit(P[Y<=2]) -2.0681 -0.92898 0.33184 0.89296 1.8375

Coefficients:


(Intercept):1 -1.38881 0.11410 -12.1718

(Intercept):2 -0.21023 0.10896 -1.9296

housingatrium -0.23210 0.13946 -1.6643


housingtower -0.52146 0.11762 -4.4336

influencelow 1.23730 0.12545 9.8631

influencemedium 0.68891 0.12302 5.6000


Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2])

Dispersion Parameter for cumulative family: 1




Para verificar se a covariavel influence e significativa podemos utilizar o teste da razaode verossimilhanca:

> pchisq(2 * (logLik(fit2) - logLik(fit1)), df = length(coef(fit2)) -

+ length(coef(fit1)), lower.tail = FALSE)

[1] 9.835079e-23

Premissas do modelo de chances proporcionais

Premissa de ordem

Uma premissa basica dos modelos de regressao ordinal e que a variavel resposta Y apresenteuma certa ordem de acordo com cada covariavel X. Assumindo que a covariavel X se


relaciona linearmente com o logaritmo da chance de Y , uma forma simples de avaliar a ordemda variavel resposta Y e verificar o grafico da media de X estratificado pelas categorias de Y .Essas medias devem apresentar uma certa ordem consistente. Se para algumas covariaveisX ′s, duas categorias adjacentes de Y nao distinguir as medias, temos evidencia de que taiscategorias de Y devem ser colapsadas.

Premissa de proporcionalidade

Juntamente com a premissa de ordem, podemos avaliar a premissa de proporcionalidadeatraves do calculo do valor esperado de X estratificado por Y , E(X|Y = j), assumindo ummodelo de chances proporcionais. Se X for uma variavel discreta, entao:

Pr(X|Y = j) =Pr(Y = j|X)Pr(X)

P (Y = j)(9)

E(X|Y = j) =∑

x

xPr(Y = j|X)Pr(X)/P (Y = j) (10)

Assumindo um modelo de chances proporcionais para Pr(Y = j|X), o valor esperadoE(X|Y = j) pode ser estimado por:

E(X|Y = j) =∑

x

xP r(Y = j|X)fx/gj (11)

em que P r(Y = j|X) sao os valores ajustados do modelo de chances proporcionais, fx e afrequencia de X = x e gj e a frequencia de Y = j. Se X for uma covariavel contınua, temos:

E(X|Y = j) =∑

x

xP r(Y = j|X)/gj (12)

A Figura abaixo ilustra o modelo de chances proporcionais para uma variavel respostaY com 4 categorias (j = 1, 2, 3, 4) e uma covariavel X. Exemplo de um modelo hipoteticocom X variando entre −50 e 50 e com β = 1 e αj = 30, 0,−10


x

Pro

babi

lidad

e

0.0

0.2

0.4

0.6

0.8

1.0

P(Y<=1) P(Y<=2) P(Y<=3)

No R podemos usar a funcao plot.xmean.ordinaly() para avaliar as premissas deordem e de proporcionalidade. Nos graficos abaixo, a linha contınua representa as mediasde X para cada categoria de Y calculada de forma empırica. Se a premissa de ordem erespeitada, sera observado uma tendencia monotonica.

As linhas pontilhadas sao os valores esperados de X|Y assumindo que a chance e propor-cional. Espera-se que as linhas contınuas e pontilhadas sejam coincidentes caso a proporci-onalidade seja verdadeira.

> par(mfrow = c(2, 2))

> plot.xmean.ordinaly(satis2 ~ housing, data = coindiv,

+ topcats = 3)


satis2

hous

ing=

apar

tmen

ts

1 2 3

0.43

0.45

0.47

n=1681satis2

hous

ing=

tow

er

1 2 30.

180.

220.

260.

30

n=1681

satis2

hous

ing=

terr

aced

1 2 30.10

0.14

0.18

0.22

n=1681


> plot.xmean.ordinaly(satis2 ~ influence, data = coindiv,

+ topcats = 2)


satis2

influ

ence

=m

ediu

m

1 2 3

0.37

0.40

n=1681

satis2

influ

ence

=lo

w

1 2 3

0.30

0.45

n=1681

> plot.xmean.ordinaly(satis2 ~ contact, data = coindiv)


satis2

cont

act=

high

1 2 3

0.54

0.55

0.56

0.57

0.58

0.59

0.60

n=1681

Resıduos

Os resıduos parciais sao os mais indicados na avaliacao de modelos de regressao para variaveisordinais. O resıduo parcial para o modelo binario e dado por:

rip = βpxip +Yi − Pi

Pi(1− Pi)(13)

em que

Pi =1

1 + exp[−(αj + x′

iβ)]

O grafico suavisado (loess) de xip contra rip fornece uma estimativa nao parametrica de comoXp se relaciona com o logaritmo da chance de Y = 1|Xp. Para o modelo ordinal e necessariocalcular o resıduo parcial para todos os Y = j, tal que:

rip = βpxip +[Yi ≥ j]− Pi

Pi(1− Pi)(14)

O grafico a ser analisado e composto por curvas suavisadas dos residuos parciais para todos osnıveis j da variavel resposta Y . Se o modelo estiver bem ajustado entao espera-se encontrar


formas e inclinacoes semelhantes para uma dada covariavel Xp para todos as categorias j.Os graficos de resıduos parcias permitem avaliar a transformacao da covariavel (linearidade)e ao mesmo tempo permite avaliar a premissa de proporcionalidade (paralelismo).

Para calcular os resıduos parciais, temos que ajustar o modelo de chances proporcionaisatraves da funcao lmr() da biblioteca Hmisc. Nao se espantem com os sinais trocados doscoeficientes comparados com os resultados da funcao vglm, pois a funcao lmr() estima omodelo para Pr(Y ≥ j) e o resultado anterior foi estimado para Pr(Y ≤ j)

> fit2a <- lrm(satis2 ~ housing + influence, x = TRUE,

+ y = TRUE, data = coindiv)

> fit2a

Logistic Regression Model

lrm(formula = satis2 ~ housing + influence, data = coindiv, x = TRUE,

y = TRUE)

Frequencies of Responses

1 2 3

567 446 668

Obs Max Deriv Model L.R. d.f. P

1681 1e-11 155.42 5 0

C Dxy Gamma Tau-a R2

0.641 0.282 0.314 0.186 0.1

Brier

0.208

Coef S.E. Wald Z P

y>=2 1.3888 0.1149 12.08 0.0000

y>=3 0.2102 0.1098 1.91 0.0555

housing=atrium 0.2321 0.1375 1.69 0.0914

housing=terraced -0.4926 0.1329 -3.71 0.0002

housing=tower 0.5215 0.1182 4.41 0.0000

influence=low -1.2373 0.1260 -9.82 0.0000

influence=medium -0.6889 0.1236 -5.57 0.0000


> resid(fit2a, "partial", pl = T)


0.0 0.4 0.8

−0.

10.

00.

10.

20.

30.

4

housing=atrium

Par

tial R

esid

ual

y>=2

y>=3

0.0 0.4 0.8

−0.

6−

0.4

−0.

20.

0

housing=terraced

Par

tial R

esid

ual

y>=2

y>=3

0.0 0.4 0.8

−0.

10.

10.

30.

5

housing=towerP

artia

l Res

idua

l y>=2

y>=3

0.0 0.4 0.8

−1.

4−

1.0

−0.

6−

0.2

influence=low

Par

tial R

esid

ual

y>=2

y>=3

0.0 0.4 0.8

−0.

8−

0.6

−0.

4−

0.2

0.0

influence=medium

Par

tial R

esid

ual

y>=2

y>=3

Modelo de razao contınua

O modelo de razao contınua baseia-se nas probabilidade condicionais em contraste com omodelo anterior que baseia-se nas probabilidade cumulativas. Seja Y = j para j = 1, . . . , ka variavel resposta. O modelo de razao contınua para Y e definido pelas seguintes k − 1regressoes:

Pr(Yi = j|Yi > j) =1

1 + exp[−(αj + x′

iβ)](15)

logit(Pr(Yi = j|Yi > j)) = ln

(Pr(Yi = j)

Pr(Yi > j)

)= α0 + x′

iβ

O modelo de razao contınua e indicado para modelar uma variavel resposta ordinalquando o indivıduo tem que passar por uma categoria para chegar ate a proxima.


Premissas do modelo de razao contınua

O modelo de razao contınua assume que o vetor de coeficientes de regressao β e o mesmoindependente de qual probabilidade condicional esta sendo estimada. Assim, assume-se umasimples razao de chances para todo Y > j, j = 1, . . . , k − 1. Se a linearidade e garantida, arazao de chances e dada por exp(βp) qualquer que seja o evento condicional Y > j.

Voltando ao exemplo, vamos estimar o modelo de risco contınuo

> fit3 <- vglm(satis2 ~ housing + influence, sratio(parallel = T),

+ coindiv)

> summary(fit3)

Call:

vglm(formula = satis2 ~ housing + influence, family = sratio(parallel = T),

data = coindiv)

Pearson Residuals:


logit(P[Y=1|Y>=1]) -1.1227 -0.70698 -5.1963e-01 1.1153 2.4157

logit(P[Y=2|Y>=2]) -2.0443 -0.78283 -3.8706e-05 1.1056 2.1587

Coefficients:


(Intercept):1 -1.30922 0.10564 -12.3935

(Intercept):2 -0.92600 0.10650 -8.6947

housingatrium -0.16282 0.12280 -1.3259


housingtower -0.45491 0.10442 -4.3565

influencelow 1.09090 0.11145 9.7883

influencemedium 0.61572 0.11031 5.5819


Names of linear predictors:

logit(P[Y=1|Y>=1]), logit(P[Y=2|Y>=2])

Dispersion Parameter for sratio family: 1





Exercıcio

This case study was taken form Harrel´s book (chapter 14). This case study was conductedby WHO in which vital signs and a large number of clinical signs and symptoms were used todevelop a predicitve model for an ordinal response. The data consist of laboratory assessmentof diagnosis and severity of illness related to pneumonia, meningitis and sepsis. The sampleconsisted of 4552 infants aged 90 days or less. The following laboratory data are used in theresponse variable:

• cerobrospinal fluid (CSF) cultire from a lumbar puncture (LP)

• blood culture (BC)

• arterial oxygen saturation (SaO2, a measure of lung dysfunction)

• chest X-ray (CXR)

The table below describe the ordinal outcome used to model the data.

Y Definition0 None of the below1 90% ≤ SaO2 < 95% or CXR+2 BC+ or CXR+ or SaO2 < 90%

Tabela 1: Ordinal outcome scale Y definition

Forty-seven clinical signs were collected for each infant (see figure below).


The signs are organized into cluster as describe below:

Task

1. Import data into R

2. Summarize/describe the Sc object


3. Y is the ordinal variable in the Sc object. Consider the covariates age, temp, rr,

hrat, waz, bul.conv, drowsy, agitated, reffort, ausc, feeding, abdominal to develop apredictive model for Y

4. Examine the assumption of ordinality of Y for each predictor by assessing how varyingY relate to the mean X, and whether the trend is monotonic. Check also the proporti-onality assumption for each covariate. The covariate temp, rr and hrat are expected tohave strong nonlinear effect. So for those predictors plot the mean absolute differencesfrom a suitable ”normal”values as an approximate solution. Use 37, 60 and 125 fortemp, rr and hrat, respectively.

5. Fit a model with all covariates in the orginal scale.

6. Check if linearity holds for all covariates by plotting the partial residuals of the modelfitted before

7. Fit a model with all covariates but transform the temp, rr and hrat covariates as before.

8. Check if there are some covariates which could be ignored to predict Y . If so, removethen and fit a reduced model.

9. Check the proportionaly assumption using the partial residuals of the fitted model.Drop the covariates of the model that do not respect this assumption.

10. Categorize the age into 3 groups using the values 6 and 70 days as cut-off points.

11. Test the interactions of age categorized and temp, hrat and rr, transformed.

12. Propose a final model and then check the partial residuals.

13. Calculate the estimated probabilities Pr(Y = j)

14. Write down a model interpretation.

R code:

load("ari.sav")

Sc$y<-Y

summary(Sc)

Sc$agecat<-cut2(Sc$age,c(7,60))

table(Sc$agecat)

options(device.ask.default=T)


plot.xmean.ordinaly(Y~age+abs(temp-37) + abs(rr-60)+

abs(hrat-125)+waz+bul.conv+drowsy+

agitated+reffort+ausc+

feeding+abdominal,data=Sc)

Documents

Modelo Linear Generalizado - Distribuição Multinomial