Valeska Andreozzi 1

Modelo Linear Generalizado

Distribuicao de Poisson

Problema 1

O objetivo desta aula e exemplificar a modelagem de dados de contagem. Vamos ilustrar como osmodelos lineares generalizados podem ser utilizados para estimar razoes de taxas.

Distribuicao de Poisson

Para dados de contagem assumimos que a variavel resposta segue uma distribuicao de Poisson comparametro µ dada pela formula:

Pr(y;µ) =µy exp(−µ)

y!Y = 0, 1, · · · ,∞ (1)

Teoricamente uma variavel aleatoria com distribuicao de Poisson pode assumir qualquer valorinteiro maior ou igual a zero. Da equacao 1 temos que a probabilidade de Y ser igual a 5 e igual a:

Pr(y;µ) =µ5 exp(−µ)

5!=

µ5 exp(−µ)

120(2)

Observe que a probabilidade de Y = 5 depende do parametro µ, que por sua vez pode dependerde covariaveis.

A distribuicao de Poisson tem como caracterısticas que E(Y ) = V ar(Y ) = µ e e geralmenteutilizada para ocorrencia de eventos raros. Da teoria estatıstica temos que a distribuicao de Poissone uma boa aproximacao da distribuicao binomial para eventos raros (numero de ensaios tende parainfinito e probabilidade sucesso tende para zero).

Considere o exemplo: Estudar a relacao entre o risco de doenca isquemica do coracao (DIC)e diversos indicadores socioeconomicos (nıvel ecologico), tendo como unidade de analise os 153bairros do Rio de Janeiro.

• Lendo o banco cardioRio.dat

> rio <- read.table("cardioRio.dat",header=T)

• Dicionario das variaveis

– bairro - nome do bairro do Rio de Janeiro

– pfave - proporcao da populacao que vive em favelas no bairro

– prede - proporcao de casas ligadas a rede publica de agua

Valeska Andreozzi 2

– pesgred - proporcao de casas ligadas a rede publica de esgotos

– pcaluga - proporcao de casas alugadas

– plixocol - proporcao de casas com coleta regular de lixo

– pesc1g - proporcao de chefes de famılia com primeiro grau completo

– palftot - proporcao da populacao que e alfabetizada

– obt3070 - numero de obitos por doenca isquemica do coracao entre 30 e 70 anos

– pop3070 - populacao entre 30 e 70 anos

– rndm2sm - proporcao de chefes de famılia com renda media ate 2 salarios mınimos

– rndm15sm - proporcao de chefes de famılia com renda media acima de 15 salarios mıni-mos

• Listando o numero de obitos e a populacao dos primeiros 10 bairros podemos observar queo numero de casos e pequeno em relacao a populacao, indicando um evento raro e que adistribuicao da variavel resposta Y pode ser a Poisson.

> rio[1:10,c("bairro","obt3070","pop3070")]

bairro obt3070 pop3070

1 SAUDE 1 1133

2 GAMBOA 7 4934

3 SANTO CRISTO 9 5737

4 CAJU 6 6935

5 CENTRO 43 26149

6 CATUMBI 13 5587

7 RIO COMPRIDO 16 18518

8 CIDADE NOVA 5 3639

9 ESTACIO 17 8273

10 SAO CRISTOVAO 30 19264

Modelando taxa

A distribuicao de Poisson tambem e muito utilizada para modelar taxas. Suponha que a taxa demortalidade de DIC por tempo de observacao seja dada por

λ =µ

l(3)

em que µ e o numero de eventos esperados (numero de casos de uma determinada doenca) e l e aquantidade total de pessoa-tempo (exemplo: tempo total de pessoas em risco de ter a doenca) emcada subgrupo de interesse. No exemplo da DIC, temos uma estimativa de risco ao inves de taxa,pois l e igual a populacao de cada bairro.

Valeska Andreozzi 3

Suponha que os dados do exemplo da DIC tenham sido obtidos atraves do acompanhamento aolongo de um ano em que, ao inves de populacao, tivessemos pessoa-tempo e quisessemos modelara taxa de DIC. Neste caso terıamos o seguinte modelo:

log(λi) = β0 + β1xi1 + β2xi2

log

(

µi

li

)

= β0 + β1xi1 + β2xi2

log(µi)− log(li) = β0 + β1xi1 + β2xi2

log(µi) = log(li) + β0 + β1xi1 + β2xi2 (4)

Tecnicamente, assumindo um processo de Poisson com intensidade λ, tem-se para a contagemY no intervalo de tempo l uma distribuicao de Poisson Yi ∼ Poisson(liλi). Consequentemente, amedia µi = liλi depende do intervalo de tempo li. Seja a dependencia da intensidade λ em relacaoas covariaveis da forma log-linear (log(λi)), entao o E(Yi) = µi e dado por

log(µi) = log(li) + β0 + β1xi1 + β2xi2 (5)

A equacao (4) ilustra a forma de modelar a taxa tendo como dados a contagem e pessoa-tempo(ou populacao). Este ultimo dado e conhecido como offset, pois nao lhe e atribuıdo nenhumcoeficiente a ser estimado. Resumindo, temos:

Yi ∼ Poisson(µi)

λi =µi

lilog(µi) = log(li) + β0 + β1xi1 + β2xi2

• Faca uma analise exploratoria das variaveis do banco.

> summary(rio)

bairro pfave prede

ABOLICAO : 1 Min. :0.00000 Min. :0.0000

ACARI : 1 1st Qu.:0.02504 1st Qu.:0.9706

AGUA SANTA : 1 Median :0.08928 Median :0.9875

ALTO DA BOA VISTA: 1 Mean :0.15992 Mean :0.9594

ANCHIETA : 1 3rd Qu.:0.21715 3rd Qu.:0.9960

ANDARAI : 1 Max. :0.90407 Max. :1.0000

(Other) :147

pesgred pcaluga plixocol

Min. :0.00191 Min. :0.0000 Min. :0.1364

Valeska Andreozzi 4

1st Qu.:0.74042 1st Qu.:0.1886 1st Qu.:0.9393

Median :0.89767 Median :0.2483 Median :0.9797

Mean :0.73820 Mean :0.2562 Mean :0.9461

3rd Qu.:0.96575 3rd Qu.:0.3131 3rd Qu.:0.9949

Max. :1.00000 Max. :0.6683 Max. :1.0000

pesc1g palftot obt3070

Min. :0.0909 Min. :0.6837 Min. : 0.00

1st Qu.:0.4114 1st Qu.:0.8873 1st Qu.: 5.00

Median :0.5018 Median :0.9208 Median : 12.00

Mean :0.5307 Mean :0.9115 Mean : 18.07

3rd Qu.:0.6864 3rd Qu.:0.9476 3rd Qu.: 24.00

Max. :0.8935 Max. :0.9772 Max. :125.00

pop3070 rndm2sm rndm15sm

Min. : 31 Min. :0.0556 Min. :0.0009

1st Qu.: 5272 1st Qu.:0.2637 1st Qu.:0.0111

Median : 10588 Median :0.4230 Median :0.0312

Mean : 15899 Mean :0.3993 Mean :0.0807

3rd Qu.: 20925 3rd Qu.:0.5220 3rd Qu.:0.0897

Max. :100804 Max. :0.7358 Max. :0.5594

> hist(rio$obt3070,main="Obitos",xlab="obitos",ylab="frequencia")

Valeska Andreozzi 5

Óbitos

óbitos

freq

uênc

ia

0 20 40 60 80 100 120

010

2030

4050

6070

Multicolinearidade

Multicolinearidade e o grau de correlacao existente entre as covariaveis. Uma correlacao forte entrecovariaveis produz grande variabilidade nas estimativas dos coeficientes de regressao. Especifica-mente, os coeficientes podem mudar drasticamente dependendo que termos estao dentro ou fora domodelo ou em que ordem eles foram introduzidos no modelo.

Uma forma de avaliar a multicolinearidade e atraves do grafico de dispersao ou da matriz decorrelacao. Uma outra alternativa e calcular o VIF (Variance Inflation factor). O VIF forneceuma medida de quanto a variancia da estimativa dos coeficientes e inflacionada comparado quandoas covariaveis nao estao linearmente dependente.

V IFp =1

1−R2p

em que R2p e um coeficiente de determinacao multipla da regressao da covariavel Xp em todas as

outras covariaveis. Suponha 3 covariaveis, X1, X2, X3. R2

1e igual ao coeficiente de determinacao

da regressao X1 ∼ X2 +X3, e assim sucessivamente.

Quando V IFp ≈ 1, isto e, R2p ≈ 0, temos que as covariaveis sao independentes e quando V IFp

e maior que 10 implica que as covariaveis estao linearmente dependente (este ponto de corte earbitrario). A raiz quadrada de V IFp pode ser interpretada como uma aproximacao de quantas

Valeska Andreozzi 6

vezes o erro padrao da covariavel Xp esta aumentado comparado com o seu erro padrao caso naohouvesse colinearidade.

O que fazer quando multicolinearidade esta presente:

1. Ignorar o problema. Quanto o objetivo da analise e predicao, os resultado devem ser adequa-dos.

2. Aumentar o tamanho da amostra, principalmente se os dados sao poucos. Isto pode reduzira correlacao entre as covariaveis.

3. Nao considerar algumas variaveis e ajustar um modelo mais simples.

4. Recodificar a covariavel ou usar uma proxy.

• Avalie a multicolinearidade entre as covariaveis atraves do diagrama de dispersao

> plot(rio[,c(2:8,11:12)])

pfave

0.0 0.8 0.0 0.5 0.2 0.8 0.1 0.6

0.0

0.6

0.0

0.8

prede

pesgred

0.0

0.8

0.0

0.5

pcaluga

plixocol

0.2

0.8

0.2

0.8

pesc1g

palftot

0.70

0.95

0.1

0.6

rndm2sm

0.0 0.6 0.0 0.8 0.2 0.8 0.70 0.95 0.0 0.4

0.0

0.4

rndm15sm

• Verifique a correlacao entre as covariaveis

> correlacao<-cor(rio[,c(2:8,11:12)])

> round(correlacao,digits=2)

Valeska Andreozzi 7

pfave prede pesgred pcaluga plixocol pesc1g palftot

pfave 1.00 -0.07 -0.11 -0.29 -0.09 -0.56 -0.65

prede -0.07 1.00 0.40 0.37 0.82 0.34 0.57

pesgred -0.11 0.40 1.00 0.50 0.53 0.42 0.49

pcaluga -0.29 0.37 0.50 1.00 0.41 0.29 0.46

plixocol -0.09 0.82 0.53 0.41 1.00 0.49 0.66

pesc1g -0.56 0.34 0.42 0.29 0.49 1.00 0.85

palftot -0.65 0.57 0.49 0.46 0.66 0.85 1.00

rndm2sm 0.58 -0.24 -0.35 -0.21 -0.41 -0.97 -0.79

rndm15sm -0.25 0.02 0.16 -0.08 0.16 0.74 0.40

rndm2sm rndm15sm

pfave 0.58 -0.25

prede -0.24 0.02

pesgred -0.35 0.16

pcaluga -0.21 -0.08

plixocol -0.41 0.16

pesc1g -0.97 0.74

palftot -0.79 0.40

rndm2sm 1.00 -0.73

rndm15sm -0.73 1.00

• Calcule o VIF para avaliar a multicolinearidade

> library(rms)

> vif(glm(obt3070~pfave+prede+pesgred+pcaluga+plixocol+

+ pesc1g+palftot+rndm2sm+rndm15sm, family=poisson, data=rio))

pfave prede pesgred pcaluga plixocol pesc1g

2.699517 2.055002 1.707209 1.587473 2.576274 38.503206

palftot rndm2sm rndm15sm

10.903771 24.841433 4.293883

> vif(glm(obt3070~pfave+prede+pesgred+pcaluga+plixocol+

+ palftot+rndm2sm+rndm15sm, family=poisson, data=rio))

pfave prede pesgred pcaluga plixocol palftot rndm2sm

2.676341 2.062273 1.725117 1.503600 2.612693 8.579149 7.503678

rndm15sm

2.841414

• Qual seria sua hipotese para relacionar obito por DIC e situacao socioeconomica?Resp.:Parece que ao retirar a variavel escolaridade pesc1g resolve o problema de colineari-dade.

Valeska Andreozzi 8

• A variavel prede (proporcao de casas ligadas a rede publica de agua) e quase toda acima de95%. Sugerimos usa-la como covariavel categorica

> rio$redecat <- ifelse(rio$prede>.95,1,0)

• Por que devemos modelar a taxa de mortalidade de DIC ao inves do numero de obitos porDIC?Resposta: Devemos modelar a taxa pois a populacao dos bairros sao diferentes, logo 10 obitosnum bairro com 10 mil habitantes e diferente de 10 obitos ocorridos num bairro com umapopulacao de mil habitantes

• E vamos ao modelo linear generalizado.

• Para covariaveis com correlacao alta, escolha somente uma delas para incluir no modelo. Suaescolha pode basear-se:

– na importancia epidemiologica das covariaveis,

– nos resultados dos modelos de regressao bivariados

– nas estimativas dos coeficientes de regressao (β′s) com maior significancia estatıstica

– no VIF

• Para modelar a taxa devemos incluir no modelo um offset igual ao logaritmo da populacao.

> rio.glm1 <- glm(obt3070~ pfave + offset(log(pop3070)),

+ data=rio, family=poisson)

• Experimente outros modelos com uma so variavel explicativa. Olhando a tabela abaixo, qualvariavel voce acrescentaria primeiro no modelo? Pare este efeito sera criado uma funcaoque resulta na analise de deviance de cada modelo com uma so covariavel comparado com omodelo nulo.

Criando uma funcao para montar esta tabela. A funcao deviancef tem como argumento onumero da coluna da covariavel que sera acrescentada no modelo referencia.

> deviancef<-function(x){

+ modelref<-glm(obt3070~ offset(log(pop3070)), data=rio, family=poisson)

+ modelo<-glm(obt3070~rio[,x]+offset(log(pop3070)), data=rio, family=poisson)

+ teste<-anova(modelref,modelo,test="Chisq")

+ resp<-teste[2,]

+ resp

+ }

> variaveis<-c(2,4:8,13)

> x<-t(sapply(variaveis,deviancef))

> resultado<-data.frame(names(rio[variaveis]),x)

> names(resultado)<-c("Variavel",attributes(x)$dimnames[[2]])

> resultado

Valeska Andreozzi 9

Variavel Resid. Df Resid. Dev Df Deviance Pr(>Chi)

1 pfave 151 263.5179 1 26.12927 3.193064e-07

2 pesgred 151 273.0289 1 16.61833 4.570718e-05

3 pcaluga 151 233.0964 1 56.55083 5.476422e-14

4 plixocol 151 286.9938 1 2.653375 0.1033303

5 pesc1g 151 281.6482 1 7.999023 0.00468026

6 palftot 151 265.4407 1 24.20655 8.653754e-07

7 redecat 151 286.6368 1 3.010408 0.08273147

A funcao seguinte resulta na analise de deviance ao se acrescentar mais uma covariavel nomodelo que possui somente a covariavel pcaluga, pois foi essa variavel que na analise anteriorreduziu mais a funcao desvio e de forma significativa.

> deviancef<-function(x){

+ modelref<- glm(obt3070~pcaluga+ offset(log(pop3070)), data=rio, family=poisson)

+ modelo<-glm(obt3070~rio[,x]+pcaluga+offset(log(pop3070)), data=rio, family=poisson)

+ teste<-anova(modelref,modelo,test="Chisq")

+ resp<-teste[2,]

+ resp

+ }

>

> variaveis<-c(2,4,6:8,13)

> x<-t(sapply(variaveis,deviancef))

> resultado<-data.frame(paste(names(rio[variaveis]),"+ pcaluga"),x)

> names(resultado)<-c("Variavel",attributes(x)$dimnames[[2]])

> resultado

Variavel Resid. Df Resid. Dev Df Deviance

1 pfave + pcaluga 150 220.6995 1 12.39685

2 pesgred + pcaluga 150 233.0948 1 0.001559777

3 plixocol + pcaluga 150 231.4248 1 1.671561

4 pesc1g + pcaluga 150 233.0363 1 0.06010507

5 palftot + pcaluga 150 230.2072 1 2.88921

6 redecat + pcaluga 150 232.9191 1 0.177295

Pr(>Chi)

1 0.0004300592

2 0.9684965

3 0.1960495

4 0.8063299

5 0.08917464

6 0.6737083

• Na tabela de analise de deviance, temos que a variavel pfave e a unica que contribui paraum aumento do deviance de forma significativa.

Valeska Andreozzi 10

> rio.glm2 <- glm(obt3070~pcaluga+pfave+ offset(log(pop3070)), data=rio, family=poisson)

> summary(rio.glm2)

Call:

glm(formula = obt3070 ~ pcaluga + pfave + offset(log(pop3070)),

family = poisson, data = rio)

Deviance Residuals:

Min 1Q Median 3Q Max

-3.8229 -0.7889 -0.0845 0.5906 3.1231

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -7.06331 0.06193 -114.053 < 2e-16 ***

pcaluga 1.35920 0.20363 6.675 2.47e-11 ***

pfave -0.44127 0.12826 -3.440 0.000581 ***

---

Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 289.65 on 152 degrees of freedom

Residual deviance: 220.70 on 150 degrees of freedom

AIC: 863.72

Number of Fisher Scoring iterations: 4

• Vamos agora interprete o modelo acima estimado.

As covariaveis pcaluga e pfave sao proporcoes que variam entre 0 e 1. Logo, devemos levarem conta esse domınio da covariavel na intrepretacao do modelo.

Para facilitar o processo de interpretacao, vamos escrever a equacao do modelo ajustado:

Yi = obitos por DIC ∼ Pois(λi)

E(Yi) = λi

V ar(Yi) = λi

λi =µi

popi

ln(λi) = −7.06331 + 1.35920pcalugai − 0.44127pfavei

Podemos dizer que para um aumento de 10% na proporcao de casas alugas, a taxa de obitopor DIC aumenta 14% (exp(β1 × 0.10) = exp(1.35920× 0.10) = 1.14).

Valeska Andreozzi 11

Ja um aumento na proporcao da populacao que vive em favelas, temos uma reducao na taxade obito por DIC de 5% (exp(β2 × 0.10) = exp(−0.44127 ∗ .10) = 0.95).

• Poderıamos continuar selecionando variaveis da forma anterior (analise de deviance) ou uti-lizar o processo stepwise.

• Para facilitar a utilizacao da funcao step(), criamos um novo objeto (rio2) que nao possuio nome do bairro e a variavel prede

> rio2 <- rio[,c(2,4:13)]

> rio.glm <-glm(obt3070~pfave + pesgred + pcaluga + plixocol + pesc1g +

+ palftot + rndm2sm + rndm15sm + redecat + offset(log(pop3070)),

+ data=rio2, family=poisson)

> rio.glm3 <- step(rio.glm, direction="both")

Start: AIC=854.26

obt3070 ~ pfave + pesgred + pcaluga + plixocol + pesc1g + palftot +

rndm2sm + rndm15sm + redecat + offset(log(pop3070))

Df Deviance AIC

- redecat 1 197.24 852.26

- palftot 1 197.41 852.43

- plixocol 1 198.69 853.70

- pesgred 1 198.69 853.71

<none> 197.24 854.26

- pfave 1 202.14 857.16

- rndm2sm 1 203.59 858.61

- pesc1g 1 205.59 860.61

- rndm15sm 1 207.05 862.07

- pcaluga 1 210.08 865.10

Step: AIC=852.26

obt3070 ~ pfave + pesgred + pcaluga + plixocol + pesc1g + palftot +

rndm2sm + rndm15sm + offset(log(pop3070))

Df Deviance AIC

- palftot 1 197.41 850.43

- pesgred 1 198.69 851.71

- plixocol 1 198.93 851.95

<none> 197.24 852.26

+ redecat 1 197.24 854.26

- pfave 1 202.15 855.17

- rndm2sm 1 203.62 856.64

- pesc1g 1 205.59 858.61

- rndm15sm 1 207.05 860.07

Valeska Andreozzi 12

- pcaluga 1 210.14 863.16

Step: AIC=850.43

obt3070 ~ pfave + pesgred + pcaluga + plixocol + pesc1g + rndm2sm +

rndm15sm + offset(log(pop3070))

Df Deviance AIC

- pesgred 1 198.71 849.73

<none> 197.41 850.43

- plixocol 1 200.10 851.12

+ palftot 1 197.24 852.26

+ redecat 1 197.41 852.43

- rndm2sm 1 203.76 854.78

- pfave 1 204.17 855.19

- pesc1g 1 206.08 857.10

- rndm15sm 1 207.80 858.82

- pcaluga 1 210.15 861.17

Step: AIC=849.73

obt3070 ~ pfave + pcaluga + plixocol + pesc1g + rndm2sm + rndm15sm +

offset(log(pop3070))

Df Deviance AIC

- plixocol 1 200.65 849.67

<none> 198.71 849.73

+ pesgred 1 197.41 850.43

+ palftot 1 198.69 851.71

+ redecat 1 198.71 851.73

- pfave 1 204.25 853.27

- rndm2sm 1 205.31 854.33

- pesc1g 1 208.79 857.81

- rndm15sm 1 210.43 859.45

- pcaluga 1 216.09 865.11

Step: AIC=849.67

obt3070 ~ pfave + pcaluga + pesc1g + rndm2sm + rndm15sm + offset(log(pop3070))

Df Deviance AIC

<none> 200.65 849.67

+ plixocol 1 198.71 849.73

+ palftot 1 200.02 851.04

+ pesgred 1 200.10 851.12

+ redecat 1 200.13 851.15

- pfave 1 207.29 854.31

Valeska Andreozzi 13

- rndm2sm 1 208.57 855.59

- pesc1g 1 210.88 857.90

- rndm15sm 1 211.11 858.13

- pcaluga 1 216.44 863.45

> summary(rio.glm3)

Call:

glm(formula = obt3070 ~ pfave + pcaluga + pesc1g + rndm2sm +

rndm15sm + offset(log(pop3070)), family = poisson, data = rio2)

Deviance Residuals:

Min 1Q Median 3Q Max

-3.0126 -0.7750 0.0102 0.6276 3.2604

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -9.3239 0.7757 -12.019 < 2e-16 ***

pfave -0.4186 0.1638 -2.555 0.01063 *

pcaluga 0.9941 0.2446 4.064 4.82e-05 ***

pesc1g 2.7908 0.8814 3.166 0.00154 **

rndm2sm 2.3664 0.8459 2.798 0.00515 **

rndm15sm -1.2060 0.3778 -3.192 0.00141 **

---

Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 289.65 on 152 degrees of freedom

Residual deviance: 200.65 on 147 degrees of freedom

AIC: 849.67

Number of Fisher Scoring iterations: 4

• Voce mudaria este modelo?

• Lembre-se da multicolinearidade? Que tal retirar uma das covariaveis de renda.

> library(car)

> scatterplotMatrix(rio[,c(2:8,11:12)])

> round(cor(rio[,c(2:8,11:12)]),digits=2)

pfave prede pesgred pcaluga plixocol pesc1g palftot

pfave 1.00 -0.07 -0.11 -0.29 -0.09 -0.56 -0.65

Valeska Andreozzi 14

prede -0.07 1.00 0.40 0.37 0.82 0.34 0.57

pesgred -0.11 0.40 1.00 0.50 0.53 0.42 0.49

pcaluga -0.29 0.37 0.50 1.00 0.41 0.29 0.46

plixocol -0.09 0.82 0.53 0.41 1.00 0.49 0.66

pesc1g -0.56 0.34 0.42 0.29 0.49 1.00 0.85

palftot -0.65 0.57 0.49 0.46 0.66 0.85 1.00

rndm2sm 0.58 -0.24 -0.35 -0.21 -0.41 -0.97 -0.79

rndm15sm -0.25 0.02 0.16 -0.08 0.16 0.74 0.40

rndm2sm rndm15sm

pfave 0.58 -0.25

prede -0.24 0.02

pesgred -0.35 0.16

pcaluga -0.21 -0.08

plixocol -0.41 0.16

pesc1g -0.97 0.74

palftot -0.79 0.40

rndm2sm 1.00 -0.73

rndm15sm -0.73 1.00

pfave

0.0 0.8 0.0 0.5 0.2 0.8 0.1 0.6

0.0

0.6

0.0

0.8 prede

pesgred

0.0

0.8

0.0

0.5 pcaluga

plixocol

0.2

0.8

0.2

0.8 pesc1g

palftot

0.70

0.95

0.1

0.6 rndm2sm

0.0 0.6 0.0 0.8 0.2 0.8 0.70 0.95 0.0 0.4

0.0

0.4rndm15sm

• Estime o modelo sem a covariavel rndm2sm. Faca as mudancas que considerar importantese compare com o modelo anterior usando analise de deviance (se os modelos se mantiveremencaixados)

Valeska Andreozzi 15

> rio.glm4 <- update(rio.glm3,~.-rndm2sm)

> anova(rio.glm4,rio.glm3,test="Chisq")

Analysis of Deviance Table

Model 1: obt3070 ~ pfave + pcaluga + pesc1g + rndm15sm + offset(log(pop3070))

Model 2: obt3070 ~ pfave + pcaluga + pesc1g + rndm2sm + rndm15sm + offset(log(pop3070))

Resid. Df Resid. Dev Df Deviance Pr(>Chi)

1 148 208.57

2 147 200.65 1 7.9203 0.004888 **

---

Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

• Escreva a hipotese nula da analise de deviance. Qual a conclusao do teste?

A Hipotese nula do teste e de que o coeficiente da variavel rndm2sm e igual a zero. O testerejeita a H0.

Se o objetivo do estudo fosse predicao, poderıamos ter optado pelo modelo rio.glm3 que in-clui covariaveis colineares. No entanto, como o objetivo e estudar associacao, vamos descartara covariavel rndm2sm.

• Vamos a analise de resıduos.

Quando o termo offset esta presente no modelo de Poisson, o R estima os preditores linearese os valores ajustados para µi ao inves de estimar para λi, ie,

preditor linear = ηi = ln popi + βxi

valores ajustados = µi = exp(ηi) = exp(ln popi + βxi)

Na verdade o que queremos sao os preditores lineares e os valores ajustados de λi. Por issotemos que fazer algumas alteracoes no resultado do R.

> res<-rstandard(rio.glm4, type = "deviance")

> pred.mu<-rio.glm4$linear.predictors #normal, poisson, gamma

> pred.lambda<-pred.mu-log(rio2$pop3070) #poisson com offset

> plot(pred.lambda,res,ylab = "Resıduo deviance padronizado",

+ xlab= "Preditor linear")

> abline(h=0)

Valeska Andreozzi 16

−7.2 −7.0 −6.8 −6.6 −6.4 −6.2

−3

−2

−1

01

23

Preditor linear

Res

íduo

dev

ianc

e pa

dron

izad

o

> fitted.mu<-rio.glm4$fitted.values

> fitted.lambda<-fitted.mu/rio2$pop3070

> pred<-2*sqrt(fitted.lambda) #poisson com offset

> plot(pred,res,ylab = "Resıduo deviance padronizado",

+ xlab ="2*sqrt(valores ajustados)")

> abline(h=0)

Valeska Andreozzi 17

0.05 0.06 0.07 0.08 0.09

−3

−2

−1

01

23

2*sqrt(valores ajustados)

Res

íduo

dev

ianc

e pa

dron

izad

o

> source("glmfunc.r")

> plotleverage(rio.glm4)

Valeska Andreozzi 18

0 50 100 150

02

46

Índice

Leve

rage

h/(

p/n)

Leverage

2

5

8

10 20

24

25

26

27

293341

51

63

128

141

144

149

151

> plotcooks(rio.glm4)

Valeska Andreozzi 19

0 50 100 150

0.0

0.1

0.2

0.3

Índice

Coo

k’s

dist

ance

Cook´s Distance

128

33

1495129

• Identifique quem sao os bairros com maior influencia no modelo

> rio[c(128,33),]

bairro pfave prede pesgred pcaluga

128 BARRA DA TIJUCA 0.01255685 0.9860039 0.6177425 0.2204779

33 TIJUCA 0.10665023 0.9632327 0.9178057 0.2384735

plixocol pesc1g palftot obt3070 pop3070 rndm2sm rndm15sm

128 0.9910596 0.8935 0.9579 15 30496 0.0767 0.5259

33 0.9581690 0.7584 0.9513 125 85542 0.1875 0.2187

redecat

128 1

33 1

• Interprete os graficos dos resıduos

Os graficos dos resıduos sugerem que o modelo apresentado nao apresenta heterocedasticidade.

• O que podemos concluir com este modelo. Ele e util?

Parece ser util.

• Observe o teste da saıda da funcao gof(). O modelo se ajusta aos dados?

> gof(rio.glm4)

Valeska Andreozzi 20

Hip. nula: Modelo e adequado

Deviance = 208.5691 com 148 Graus de Liberdade

P-valor 0.00077

Contudo parece que o modelo nao se ajusta bem aos dados. Novamente aqui ha que levar emconta o objetivo do modelo (predicao ou descricao)

• Acrescente a seguinte covariavel ao modelo: proporcao da populacao entre 60 e 70 anos

> rio$pidosos <- scan("ppidosos.dat")

• Ajuste outros modelos agora com esta covariavel indicadora da estrutura etaria da populacao

> rio.glm5<-glm(obt3070 ~ pfave + pcaluga + pesc1g + rndm15sm +

+ pidosos + offset(log(pop3070)),

+ family = poisson, data = rio)

> summary(rio.glm5)

Call:

glm(formula = obt3070 ~ pfave + pcaluga + pesc1g + rndm15sm +

pidosos + offset(log(pop3070)), family = poisson, data = rio)

Deviance Residuals:

Min 1Q Median 3Q Max

-3.3120 -0.6491 -0.0766 0.6003 2.9552

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -7.4716 0.1449 -51.547 < 2e-16 ***

pfave -0.4743 0.1646 -2.882 0.00396 **

pcaluga 0.5802 0.2614 2.219 0.02645 *

pesc1g -0.1931 0.3224 -0.599 0.54920

rndm15sm -1.0823 0.3912 -2.767 0.00566 **

pidosos 5.2057 0.8860 5.876 4.21e-09 ***

---

Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 289.65 on 152 degrees of freedom

Residual deviance: 173.40 on 147 degrees of freedom

AIC: 822.42

Number of Fisher Scoring iterations: 4

Valeska Andreozzi 21

Ao incluirmos a estrutura etaria da populacao temos algumas alteracoes nos efeitos das outrascovariaveis presentes no modelo, especialmente na variavel pcaluga e pesc1g, apesar destaultima nao apresentar efeito estatisticamente significativo.