Modelos Lineares Generalizados - Verificação do Ajuste do Modelo

Estatística de Pearson GeneralizadaTestes de Hipóteses usando a função Deviance

Modelos Lineares Generalizados -Verificação do Ajuste do Modelo

Erica Castilho Rodrigues

21 de Junho de 2013

1


Estatística de Pearson Generalizada

Testes de Hipóteses usando a função Deviance

2


Estatística de Pearson Generalizada

3


Estatística de Pearson GeneralizadaI Uma outra medida usada para verificar o ajuste do modelo.I Essa estatística é dada por

X 2p =

n∑

i=1

(yi − µi)2

Var(Yi)

onde Var(Yi) é a função de variância estimada sob omodelo que está sendo ajustado aos dados.

4


I Para o Poisson e Binomial a estatística fica

X 2p =

n∑

i=1

(oi − ei)2

ei

que é a Estatística Qui-Quadrado usual.I Essa estatistica tem a seguinte distribuição assintótica

X 2p ∼ χ2

n−p

ondeI n é o tamanho da amostra;I p é o número de parâmetros do modelo.

5


I A Deviance é mais usada do que a Estatística de PearsonGeneralizada.

I Isso acontece porque para a Deviance temos que:I seu valor sempre dimui quando acrescentamos variáveis

no modelo;I o mesmo não é verdade para a Estatística de Pearson.

6


Testes de Hipóteses usando a função Deviance

7


I Vimos que podemos fazer testes sobre o vetor β utilizandoa distribuição assintótica

b ∼ N(β, I(β−1) .

I Uma alternativa:I comparar o ajuste de dois modelos;I o modelo com a variável e o modelo sem a variável.

I Um modelo deve estar contido no outro.I A diferença deve ser apenas a variável incluída/retirada.I A distribuição de probabilidade deve ser a mesma.I A função de ligação deve ser a mesma.

8


I Vamos chamar o modelo mais simples (menos variáveis)de M0.

I O modelo mais complexo (mais variáveis) será M1.I Para o modelo M0 temos a hipótese nula de que

H0 : β = β0 =

[

β1...βq

]

.

I Para o modelo M1 temos a hipótese alternativa

H1 : β = β0 =

[

β1...βp

]

.

I Observe que q < p < n.

9


I Podemos testarH0 vs H1

usando a diferença das Deviances dos dois modelos

∆D = D0−D1 = 2 [l(bmax , y)− l(b0, y)]−2 [l(bmax , y)− l(b1, y)] .

I Se os modelos estão bem ajustados temos que

D0 ∼ χ2(n−p) D0 ∼ χ2

(n−q) .

I Portanto

∆D = D0 − D1 ∼ χ2(n−q)−(n−p) ou seja χ2

p−q .

10


I Hipóteses a serem testadas:I H0: a diferença entre M0 e M1 não é significativa;I H1: a diferença enrte os modelos é significativa.

I Se ∆D não é um valor atípico na distribuição χ2p−q:

I podemos aceitar H0 permanecer com o modelo maissimples;

I a diferença de ajuste entre os modelos não é significativa.

I H0 é rejeitada para valores grandes ou pequenos de ∆D?Grandes.

I Como fica a região crítica?I se ∆D < χ2

c não rejeitamos H0 permanecemos com omodelo M0;

I se ∆D > χ2c rejeitamos H0 e ficamos com o modelo M1.

11


I A aproximação assintótica da distribuição de ∆D é melhordo que de D.

I Se temos um parâmetro de ruído para estimar,I nem sempre a Deviance poderá ser obtida diretamente dos

dados;I precisa ainda do parâmetro de ruído.

I Vimos no caso Normal, por exemplo, que

D =

∑

i(yi − yi)2

σ2

precisamos ainda estimar σ2.I Vejamos como isso é feito no exemplo a seguir.

12


ExemploI Considere o modelo linear normal

E(Yi) = µi = xTi β .

I Já vimos que a Deviance desse modelo é dada por

D =

∑

i(yi − yi)2

σ2

I Vamos usar a seguinte notação:I yi(0) é o valor ajustado pelo modelo M0;I yi(1) é o valor ajustado pelo modelo M1.

13


ExemploI A Deviance do modelo M0 (tem q parâmetros) fica

D0 =

∑

i(yi − yi(0))2

σ2

e do modelo M1 (que tem p parâmetros)

D1 =

∑

i(yi − yi(1))2

σ2 .

I Temos ainda que

D0 ∼ χ2n−q D1 ∼ χ2

n−p ∆D ∼ χ2p−q .

14


I Para não termos que encontrar σ2 vamos usar a razão

F =∆D/(p − q)D1/(n − p)

∼ Fp−q,n−p .

I Dessa maneira, F fica

F =∆(

∑

i(yi − yi(0))2 −∑

i(yi − yi(1))2)/(p − q)(∑

i(yi − yi(1))2)/(n − p)∼ Fp−q,n−p.

I Como o σ2 é cancelado nessa razão, torna-sedesncessário estimá-lo.

I Rejeitamos H0 quando F é grande.

15


ExemploI A tabela a seguir mostra os dados do peso e a idade de

gestação de bebês em um hospital.

16


Exemplo (continuação)I A figura a seguir mostra o gráfico de dispersão entre as

duas variáveis.

17


Exemplo (continuação)I Os bebês estão divididos em dois grupos:

I masculino e feminino.

I Como podemos escrever o modelo com essas duasvariáveis?

I A variável sexo entra como Dummy.I O modelo sem interação fica

Yi = β0 + β1Xi + β2Zi + εi εi ∼iid N(0, σ2)

ondeI Yi é o peso do bebê;I Xi é idade de gestaçãoI Zi é uma indicadora que representa sexo (1 - masculino, 0 -

feminino).

18


Exemplo (continuação)I Queremos verificar a necessidade de incluir o termo de

interação.I O modelo com interação é dado por

Yi = β0 + β1Xi + β2Zi + εi + β3XiZi εi ∼iid N(0, σ2)

I Vamos denotar porI M0: modelo sem interação;I M1: o modelo com interação.

I Queremos verificar se o ganho de ajuste de M1 em relaçãoa M0 é significativo.

19


Exemplo (continuação)I A Soma dos Quadrados dos Resíduos está relacionada

com a Deviance da segunte maneira

SQE =∑

i

(yi − yi)2 = σ2D.

I Para os modelos temos que

SQE0 = 658770.8 SQE1 = 652424.5

ou seja

D0 =658770.8

σ2 D1 =652424.5

σ2 .

20


Exemplo (continuação)I Temos que n = 24 logo

F =∆(

∑

i(yi − yi(0))2 −∑

i(yi − yi(1))2)/(p − q)(∑

i(yi − yi(1))2)/(n − p)

=(SQE0 − SQE1)/(p − q))

SQE1/(n − p)=

(658770.8 − 652424.5)/(4 − 3)652424.5/(24 − 4)

= 0,19 .

21


Exemplo (continuação)

I Devemos comparar esse valor com a F1,20.I Rejeitamos H0, quando F é grande.I Fixando α = 0,05, o valor crítico dessa distribuição é dado

porFc = 4.35 poisP(F1,20 > 4,35) = 0,05 .

I A região crítica é dada porI se F < Fc , não rejeitamos H0I se F > Fc , rejeitamos H0 .

I Conclusão do teste:I Fobs = 0, 19 < 4, 35 não rejeitamos H0;I não é necessário incluir termo de interação no modelo;I conclusão: o efeito da idade no peso é o mesmo para

meninos e meninas.

22


ExemploI Um pesquisador quer verificar qual a dose ideal de

inseticida para matar insetos.I Diferentes doses são usadas para grupos de uma mesma

espécie.I Vamos usar a seguinte notação:

I di : dose do inseticida;I mi : número de insetos que receberam a dose;I yi : número de insetos mortos dentre os mi que receberam

o inseticida;I pi : proporção de insetos mortos.

23


Exemplo (continuação)I A tabela a seguir mostra os dados coletados

24


Exemplo (continuação)I O pesquisador deseja determinar quais as doses tais que

I 50% dos insetos são mortos (LD50);I 90% dos insetos são mortos (LD90).

I Podem usar esse dado para aplicação em campo.

25


Exemplo (continuação)I A figura a seguir mostra o gráfico dispersão entre:

I doses de inseticida (di ) e proporção de insetos mortos (pi ).

26


Exemplo (continuação)I O gráfico tem um aspecto sigmoidal.I Esse formato pode nos guiar na escolha da função de

ligação.I Esse tipo de ensaio é chamado de dose-resposta.I Dois aspectos devem ser considerados:

I a dose da droga (inseticida, fungicida, herbicida,medicamento);

I o indivíduo que recebe a droga (inseto, planta, fungo,paciente).

I A reposta do indivíduo é binária:I responde (1) ou não responde (0) ao tratamento.

27


Exemplo (continuação)I A resposta dependerá do nível da dosagem aplicada.I Cada indivíduo tem um nível a partir do qual responde ao

tratamento.I Esse valor é chamdo de tolerância do indivíduo.I Essa tolerância varia de um indivíduo para o outro dentro

da população.I Portanto é uma variável aleatória e vamos denotá-la por U.

28


Exemplo (continuação)I A figura seguir mostra exemplos de distribuição da

tolerância.

29


Exemplo (continuação)I Vamos denotar por f (u) a função de densidade da

tolerância.I Seja d a dose ministrada à toda população.I Quais indivíduos responderão à droga?I Aqueles tais que

U < d .

I A probabilidade de um indivíduo escolhido ao acasoresponda ao tratamento é

π(d) = P(U < d) =∫ d

−∞

f (u)du .

30


Exemplo (continuação)I Para valores pequenos de d quanto deve valer π(d)?

π(d) ≈ 0 .

I Para valores grandes de d quanto deve valer π(d)?

π(d) ≈ 1 .

I π é uma função crescente ou decrescente de d?I Crescente, quanto maior a dose maior a probabilidade de

resposta.

31


Exemplo (continuação)I No exemplo dos insetos queremos encontrar um modelo

razoável de como π(d) varia com d .I E a partir disso encontrar os valores de doses tais que

I 50% dos indivíduos respondem à droga (LD50);I 90% dos indivíduos respondem à droga (LD90).

I Seja Yi a variável aleatória que denota o número deinsetos mortos.

I Seja πi a probailidade de um inseto do i-ésimo grupomorrer.

I Qual a distribuição de Yi?

Yi ∼ Bin(πi ,mi) .

32


Exemplo (continuação)I Vamos usar a função de ligação canônica.I Qual ligação é essa? Logística.I Isso significa que:

π =1

1 + eηiou log

(

πi

1 − πi

)

= ηi .

I Vamos ajustar o seguinte modelo

Yi ∼ Bin(πi ,mi) log(

πi

1 − πi

)

= β0 + β1di .

33


Exemplo (continuação)I O script usado para ajustar o modelo foi o seguinte:

x=c(0,2.6,3.8,5.1,7.7,10.2)m=c(49,50,48,46,49,50)y=c(0,6,16,24,42,44)

dados=data.frame(x=x,y=y,m=m)

modelo=glm(cbind(y,m-y)~x, family="binomial",data=dados)

I Precisamos criar dois vetores:I um com o número de sucesos e outro com número de

fracassos.

cbind(y,m-y)

34


Exemplo (continuação)I O resumo do ajuste encontra-se a seguir.

> summary(modelo)

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

(Intercept) -3.22566 0.36992 -8.720 <2e-16 ***x 0.60513 0.06781 8.923 <2e-16 ***

Null deviance: 163.745 on 5 degrees of freedomResidual deviance: 10.258 on 4 degrees of freedomAIC: 33.479

35


Exemplo (continuação)I Qual modelo estimado?

log(

πi

1 − πi

)

= −3.22 + 0.60(di) .

I Qual interpretação do β1?I Vamos tirar a exponencial dos dois lados

(

πi

1 − πi

)

= e−3,22+0,60(di) = e−3,22∗ e0,60(di)

I O que acontece se aumentarmos di em uma unidade(

πi

1 − πi

)

= e−3,22∗ e0,60(di+1) = e−3,22

∗ e0,60(di)e0,60

a razão de chance fica multiplicada por e0,6 = 1,82.36


Exemplo (continuação)I Isso equivale a aumentar 82%.I O termo

(

πi

1 − πi

)

é denominado razão de chances (odds ratio) e mede oquanto o sucesso é mais provável que o fracasso.

I Exemplo se a razão é 5, significa que a probabilidade desucesso é 5 vezes maior que a probabilidade de fracasso.

I Conclusão: para cada aumento em uma unidade da dose,a razão de chances é multiplicada por eβ1 que nesse casoequivale a aumentar 82%.

I Esse interpretação só é possível na ligação canônica.

37


Exemplo (continuação)I Vamos encontrar agora os valores estimados das doses

letais.I Lembre que:

I LD50 dose tal que 50% dos insetos são mortos;I LD90 dose tal que 90% dos insetos são mortos.

I Vimos que o modelo estimado é dado por

log(

πi

1 − πi

)

= −3,22 + 0,60(di) .

I Vamos isolar di

di =

(

log(

πi

1 − πi

)

+ 3.22)

/0,60

38


Exemplo (continuação)I LD50 corresponde ao valor de di tal que p = 50%;I LD90 corresponde ao valor de di tal que p = 90%.I Portanto

LD50 =

(

log(

0,51 − 0,50

)

+ 3.22)

/0,60 = 5,37

LD90 =

(

log(

0,91 − 0,90

)

+ 3.22)

/0,60 = 9,03

39


Exemplo (continuação)I A figura a seguir mostra o gráfico de dispersão dos dados

com a curva ajustada sobreposta.

0 2 4 6 8 10

0.0

0.2

0.4

0.6

0.8

Dose

Pro

porç

ão

40


Exemplo (continuação)I Como podemos verificar se o modelo está bem ajustado?

Deviance.I O seguinte comando retorna a Deviance do modelo

> modelo$deviance[1] 10.25832

I Com qual distribuição de referência devemos com paraesse valor?

I Com uma distribuição χ2n−p, no nosso caso χ2

4.

41


Exemplo (continuação)I Rejeitamos H0 para valores altos da deviance.I Portatno o p-valor é dado por

P(χ24 ≥ 10.25832)

1-pchisq(10.25,4)[1] 0.03642058

I Conclusão: rejeitamos H0 e concluímos que o modelo nãoestá bem ajustado.

42


Exemplo (continuação)I Queremos agora verificar se, de fato, a dose é significativa

para explicar a resposta.I Isso equivale a comparar os modelos:

I M0: ηi = β0 (modelo nulo, só com intercepto);I M1: ηi = β0 + β1di .

I Como n = 6 os graus de liberdade dos modelos são:I M0: n − p = 6 − 1 = 5;I M1: n − p = 6 − 2 = 4.

43


Exemplo (continuação)I A tabela a seguir mostra a Deviance e os graus de

liberdade para cada um dos modelos:

Modelo Graus de Liberada Devianceηi = β0 5 163,74

ηi = β0 + β1di 4 10,26

I A diferença entre as Deviances é dada por

∆D = 163,74 − 10,26 = 153,48 .

I Sabemos queδD ∼ χ2

1 .

I Rejeitamos H0 para valores grande ou pequenos de ∆D?Grandes.

I A região crítica é do tipoI ∆D < χ2

c ⇒ não rejeitamos H0 e ficamos com o modelo M0;I ∆D > χ2

c ⇒ rejeitamos H0 e ficamos com o modelo M1.44


Exemplo (continuação)I O valor crítico da χ2

1 é 3,84, pois

P(χ21 > 3,84) = 0,05 .

0 2 4 6 8 10

0.0

0.2

0.4

0.6

0.8

1.0

1.2

x

dchi

sq(x

, 1)

Região de Rejeitção

45


Exemplo (continuação)

I Conclusão:I ∆D = 153, 48 > 3.84 ⇒ rejeitamos H0;I isso singifica que a variável explicativa deve entrar no

modelo;I o ganho ao acrescentar essa variável é expressivo.

I Poderíamos testar a inclusão de mais variáveis no modelo.

46

Documents

Modelos Lineares Generalizados - Verificação do Ajuste do Modelo