Upload
others
View
8
Download
0
Embed Size (px)
Citation preview
1
CE225 – Modelos Lineares Generalizados
2
Objetivo da disciplina
• Apresentar ao aluno a teoria e aplicações dos Modelos Lineares Generalizados,
propostos originalmente em Nelder e Wedderburn (1972), que configuram extensões dos
modelos lineares clássicos (com erros normalmente distribuídos) e que permitem
analisar a relação funcional entre um conjunto de variáveis independentes e uma variável
aleatória dependente com distribuição pertencente à família exponencial de
distribuições.
• A família exponencial contempla, dentre outras, as distribuições normal, exponencial,
gama, normal inversa, Poisson, binomial e binomial negativa.
• Na sequência são descritos alguns dos problemas que serão analisados, ao longo do
semestre, usando elementos de Modelos Lineares Generalizados.
3
Exemplos de motivação
4
Exemplo 1 – Análise da resistência de uma nova fibra sintética usada na produção de
camisas. Sabendo-se que a resistência da fibra é afetada pela quantidade de algodão
utilizada, e que a quantidade de algodão no produto final, de acordo com as características
desejadas, deve estar no intervalo de 10 a 40%, um experimento é delineado com cinco
réplicas (amostras de tecidos) para cinco diferentes especificações referentes à quantidade de
algodão.
Variável resposta: Resistência da fibra (em libras/pol2).
Variável explicativa: Porcentagem de algodão no tecido, com cinco níveis: 15, 20, 25, 30 e
35%.
5
Dados
Quadro 1 – Dados de resistência (em libras/pol2) para o experimento
de fibra sintética.
Amostra de tecido Porcentagem
de algodão 1 2 3 4 5
15 7 7 15 11 9
20 12 17 12 18 18
25 14 18 18 19 19
30 19 25 22 19 23
35 7 10 11 15 11
Objetivos:
• Avaliar o efeito da porcentagem de algodão na resistência da fibra sintética;
• Identificar a porcentagem ideal de algodão de forma a se obter máxima resistência.
6
15 20 25 30 35
10
15
20
25
Porcentagem de algodão
Res
istê
ncia
da
fibra
Figura 1 – Gráfico de dispersão para as resistências das fibras sob cinco porcentagens distintas de
algodão.
7
Exemplo 2 – Amostras de 20 insetos (Heliothis virescens - praga do algodão) foram expostas
a doses crescentes do cipermetrina, dois dias depois da emergência da pupa (Collet, 2002).
Após 72h, foram contados os insetos mortos.
Variável resposta: Número de insetos mortos.
Variáveis explicativas
• Dose de cipermetrina, com níveis 1, 2, 4, 8, 16, 32 u.m.:
• Sexo (macho ou fêmea).
8
Quadro 2 – Números de insetos mortos
em amostras de 20 insetos machos e fêmeas
submetidos a doses crescentes de
cipermetrina.
Nº insetos mortos Dose Log(Dose)
Machos Fêmeas
1,0 0 1 0
2,0 1 4 2
4,0 2 9 6
8,0 3 13 10
16,0 4 18 12
32,0 5 20 16
9
Objetivos:
• Propor um modelo que descreva o aumento na mortalidade dos insetos segundo a dose de
aplicada de cipermetrina;
• Comparar as curvas de mortalidade para insetos machos e fêmeas;
• Estimar doses letais, ou seja, doses efetivas (mortais) para uma proporção p de insetos.
10
0 1 2 3 4 5
0.0
0.2
0.4
0.6
0.8
1.0
log(dose)
Pro
porç
ão d
e in
seto
s m
orto
s
Machos
Fêmeas
Figura 2 – Gráfico da proporção de insetos mortos segundo o sexo e a dose de inseticida.
11
Exemplo 3 - Uma população de mulheres indígenas vivendo numa região próxima a Phoenix,
Arizona, foi testada para diabetes de acordo com o critério estabelecido pela Organização Mundial
de Saúde. Os dados foram coletados pelo Instituto Nacional de Diabetes e Doenças Digestivas e
Renais dos EUA. São considerados os dados referentes aos 532 registros completos.
Variáveis explicativas
npreg - número de gestações;
gli - concentração de glicose no plasma no teste de tolerância à glicose oral.
bp - pressão arterial diastólica (mm Hg).
skin - espessura da prega tricipital (mm).
bmi - índice de massa corporal (peso/altura2).
ped - função pedigree diabetes.
age - idade em anos.
Variável resposta
type – diagnóstico de diabetes de acordo com o teste de Glicemia em Jejum.
12
Notas:
O Teste Oral de Tolerância à Glicose (também conhecido como Curva Glicêmica) é feito da
seguinte maneira: a pessoa com suspeita de diabetes ingere 75g de glicose diluída em água. Após duas
horas de espera, é feita a coleta de sangue para medir a taxa de glicose. Se o resultado for igual ou
superior a 200mg/dl (miligramas por decilitro), o indivíduo é considerado portador de diabetes. Se a
glicemia estiver entre 140mg/dl e 199mg/dl, então o diagnóstico é de pré-diabetes.
13
Dados - Seis primeiras linhas da base:
npreg glu bp skin bmi ped age type
1 5 86 68 28 30.2 0.364 24 No
2 7 195 70 33 25.1 0.163 55 Yes
3 5 77 82 41 35.8 0.156 35 No
4 0 165 76 43 47.9 0.259 26 No
5 0 107 60 25 26.4 0.133 23 No
6 5 97 76 27 35.6 0.378 52 Yes
Objetivos:
• Identificar fatores associados à incidência de diabetes;
• Estabelecer um modelo preditivo para o diagnóstico de diabetes.
14
Não Sim
0
2
4
6
8
10
12
14
Diagnóstico
Gra
vide
z
Não Sim
60
80
100
120
140
160
180
200
DiagnósticoG
licos
e
Não Sim
40
50
60
70
80
90
100
110
Diagnóstico
Pre
ssão
art
eria
l (m
mH
g)
Não Sim
20
40
60
80
100
Diagnóstico
Esp
essu
ra p
rega
tri
c.
Não Sim
20
25
30
35
40
45
Diagnóstico
IMC
(kg
/m2)
Não Sim
0.0
0.5
1.0
1.5
2.0
Diagnóstico
Ped
igre
e
Não Sim
20
30
40
50
60
Diagnóstico
Idad
e (a
mos
)
Não Sim
Diagnóstico
Núm
ero
de p
acie
ntes
0
20
40
60
80
100
120
Figura 3 – Distribuição das variáveis explicativas segundo o diagnóstico de diabetes.
15
Tabela 1 – Médias, desvios padrões e estatística do teste t
(comparação de duas médias), para amostras independentes,
para as variáveis explicativas segundo o diagnóstico
Diagnóstico Variável
Não Sim
Estatística
t
Número de gestações 2,9 (2,8) 4,8 (4,0) -3,56
Glicose (oral) 113 (26) 145 (30) -7,38
Pressão diastólica 69 (11) 74 (11) -2,95
Espessura – prega tricipital 27 (11) 33 (12) -3,39
Pedigree 0,41 (0,27) 0,55 (0,36) -4,51
Idade 29 (10) 37 (11) -2,70
16
npreg
60 100 140 180 20 40 60 80 100 0.0 0.5 1.0 1.5 2.0
02468101214
6080
100120140160180200
glu
bp405060708090100110
20406080
100
skin
bmi202530354045
0.00.51.01.52.0
ped
0 2 4 6 8 12 40 60 80 100 20 30 40 20 30 40 50 60
20
3040
5060
age
Figura 4 – Matriz de gráficos de dispersão para as variáveis explicativas.
17
Exemplo 4 – Na sequência, são apresentadas as cinco primeiras linhas de um banco de dados
referente a um estudo prospectivo com 100 indivíduos de pelo menos 65 anos de idade em boas
condições físicas. O objetivo do estudo é tentar relacionar o número de quedas num período de
seis meses com as seguintes variáveis explicativas, descritas na ordem em que aparecem na base:
Variáveis explicativas:
• Intervenção – Fator com níveis ‘Educ’: educação somente; ‘Educ+Exerc’: educação e exercícios
físicos;
• Sexo – Fator com níveis ‘Fem’: feminino; ‘Masc’: masculino;
• Balanço – escore de equilíbrio do indivíduo, numa escala de 0 a 100 (quanto maior o escore,
maior o equilíbrio;
• Força – escore de força do indivíduo, numa escala de 0 a 100 (quanto maior o escore, maior a
força).
Variável resposta:
• Quedas – número de quedas no período;
18
Objetivos:
• Avaliar o efeito da intervenção na prevenção das quedas;
• Identificar características dos indivíduos associadas a um maior número de quedas.
19
Dados: 10 primeiras linhas da base:
> head(geriatra,10)
quedas intervenção sexo balanço força
1 1 Educ+Exerc Fem 45 70
2 1 Educ+Exerc Fem 62 66
3 2 Educ+Exerc Masc 43 64
4 0 Educ+Exerc Masc 76 48
5 2 Educ+Exerc Fem 51 72
6 1 Educ+Exerc Masc 73 39
7 0 Educ+Exerc Masc 40 54
8 0 Educ+Exerc Fem 66 37
9 2 Educ+Exerc Masc 80 81
10 2 Educ+Exerc Masc 56 60
20
• Análise descritiva (univariada)
> summary(geriatra)
quedas intervenção sexo balan ço força
Min. : 0.00 Educ :50 Fem :47 Min. : 13.00 Min. :18.00
1st Qu.: 1.00 Educ+Exerc:50 Masc:53 1st Qu.: 39.00 1st Qu.:52.00
Median : 3.00 Median : 51.50 Median :60.00
Mean : 3.04 Mean : 52.83 Mean :60.78
3rd Qu.: 4.00 3rd Qu.: 66.25 3rd Qu.:70.25
Max. :11.00 Max. : 98.00 Max. :90.00
21
• Análise descritiva (bivariada):
o Número de quedas vs intervenção;
> with(geriatra,describeBy(quedas, intervenção, mat = TRUE,digits=2))
item group1 vars n mean sd median trimmed mad min max range skew kurtosis se
11 1 Educ 1 50 4.52 2.40 4 4.25 1.48 1 11 10 0.89 -0.07 0.34
12 2 Educ+Exerc 1 50 1.56 1.33 1 1.43 1.48 0 5 5 0.62 -0.52 0.19
o Número de quedas vs sexo;
> with(geriatra,describeBy(quedas, sexo, mat = TRUE ,digits=2))
item group1 vars n mean sd median trimmed ma d min max range skew kurtosis se
11 1 Fem 1 47 3.47 2.49 3 3.21 1.4 8 0 11 11 0.98 0.70 0.36
12 2 Masc 1 53 2.66 2.34 2 2.35 1.4 8 0 10 10 1.18 1.08 0.32
22
o Número de quedas vs nível de balanço;
> with(geriatra,describeBy(quedas, cut(balanço,4), mat = TRUE,digits=2))
item group1 vars n mean sd median trimme d mad min max range skew kurtosis se
11 1 (12.9,34.2] 1 19 2.32 1.25 2 2.3 5 1.48 0 4 4 -0.09 -1.34 0.29
12 2 (34.2,55.5] 1 39 2.95 2.70 2 2.5 8 1.48 0 11 11 1.36 1.35 0.43
13 3 (55.5,76.8] 1 31 3.32 2.60 3 3.1 6 2.97 0 9 9 0.43 -0.92 0.47
14 4 (76.8,98.1] 1 11 3.82 2.44 3 3.6 7 1.48 0 9 9 0.68 -0.31 0.74
o Número de quedas vs nível de força;
> with(geriatra,describeBy(quedas, cut(força,4), ma t = TRUE,digits=2))
item group1 vars n mean sd median trimmed mad min max range skew kurtosis se
11 1 (17.9,36] 1 5 2.80 1.30 3 2.80 1.48 1 4 3 -0.26 -1.96 0.58
12 2 (36,54] 1 25 2.32 2.19 2 2.10 2.97 0 7 7 0.61 -0.74 0.44
13 3 (54,72] 1 50 3.30 2.55 3 2.92 1.48 0 10 10 1.10 0.25 0.36
14 4 (72,90.1] 1 20 3.35 2.60 3 3.06 1.48 0 11 11 1.10 1.45 0.58
23
0 1 2 3 4 5
-4-2
02
46
Fitted values
Res
idua
ls
Residuals vs Fitted
5293
67
-2 -1 0 1 2
-10
12
3
Theoretical Quantiles
Sta
ndar
dize
d re
sidu
als
Normal Q-Q
52
93
67
0 1 2 3 4 5
0.0
0.5
1.0
1.5
Fitted values
Sta
ndar
dize
d re
sidu
als
Scale-Location52
93
67
0.00 0.05 0.10 0.15-2
-10
12
3
Leverage
Sta
ndar
dize
d re
sidu
als
Cook's distance
0.5
Residuals vs Leverage
5293
42
Figura 5 – Gráficos para os resíduos de um modelo de regressão linear múltipla ajustado aos dados de
quedas de idosos.
24
Exemplo 5 – Dados de 4624 apólices de seguros de automóveis que registraram sinistros no
período de um ano, entre 2004 e 2005.
Variáveis explicativas:
• Valor: Valor do veículo (x10.000 dólares);
• Tipo: tipo de veículo, com níveis: o BUS: ônibus; o CONVT: conversível; o COUPE; o HBACK: hathback; o HDTOP: hardtop; o MCARA: trailer motorizado; o MIBUS: mini-ônibus; o RDSTR: roadster; o SEDAN; o STNWG: station wagon; o TRUCK: caminhão; o UTE: utilitário.
• Idade: idade do veículo, com níveis 1 (veículos mais novos), 2, 3 e 4;
• Sexo: sexo do motorista, com níveis: M (masculino) e F (Feminino);
• Area: área de residência do motorista: A, B, C, D, E e F;
• Idademot: idade do motorista, com níveis: 1 (mais novos), 2, 3, 4, 5 e 6. Variável resposta:
• Quantia: valor (somado) dos sinistros apresentados no período (variável resposta).
25
Objetivos:
• Identificar fatores associados a sinistros mais caros;
• Estabelecer um modelo preditivo, que sirva de base para estabelecer a tabela de preços para
períodos futuros.
26
Dados
• Dez primeiras linhas da base:
> head(dados2,10)
Valorcar Quantia Tipo Idade Sexo Area Idadem ot
15 1.66 669.5100 SEDAN 3 M B 6
17 1.51 806.6100 SEDAN 3 F F 4
18 0.76 401.8055 HBACK 3 M C 4
41 1.89 1811.7100 STNWG 3 M F 2
65 4.06 5434.4400 STNWG 2 M F 3
66 1.39 865.7900 HBACK 3 F A 4
96 2.66 1105.7700 STNWG 1 F F 5
99 0.50 200.0000 HBACK 4 F A 5
116 1.16 739.2300 STNWG 4 F B 2
125 3.56 3230.6000 MCARA 3 M F 4
27
• Análise descritiva (univariada):
> summary(dados2)
Valorcar Quantia Tipo Idade Sexo Area Idademot
Min. : 0.000 Min. : 200.0 SEDAN :1476 1: 825 F:2648 A:1085 1: 496
1st Qu.: 1.100 1st Qu.: 353.8 HBACK :1264 2:1259 M:1976 B: 965 2: 932
Median : 1.570 Median : 761.6 STNWG :1173 3:1362 C:1412 3:1113
Mean : 1.859 Mean : 2014.4 UTE : 260 4:1178 D: 496 4:1104
3rd Qu.: 2.310 3rd Qu.: 2091.4 HDTOP : 130 E: 386 5: 614
Max. :13.900 Max. :55922.1 TRUCK : 120 F: 280 6: 365
(Other): 201
> sum(dados2$Quantia>15000) ### Numero de apólices que geraram mais de
$15.000 em sinistros.
[1] 65
28
Todas as apólices com sinistros
Valor dos sinistros (x$10.000)
Fre
quên
cia
0 5 10 15 20 25 30 35 40 45 50 55 60
0
200
400
600
800
1000
1200
1400
1600
1800
Apólices com sinistros inferiores a $15.000
Valor dos sinistros (x$10.000)
Fre
quên
cia
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0
200
400
600
800
1000
1200
1400
1600
1800
Figura 6 – Distribuição dos valores dos sinistros gerados pelos segurados.
29
• Análise descritiva (bivariada) – Valores dos sinistros vs tipo de veículo:
> with(dados2,describeBy(Quantia, Tipo, mat = TRUE,digits=2))
item group1 vars n mean sd median trimmed mad min max range skew kurtosis se
11 1 BUS 1 9 1484.79 1483.53 876.48 1484.79 613.99 371.82 4790.84 4419.02 1.13 -0.11 494.51
12 2 CONVT 1 3 2296.27 3319.80 530.00 2296.27 440.33 233.00 6125.81 5892.81 0.38 -2.33 1916.69
13 3 COUPE 1 68 2760.64 4197.29 1171.92 1819.55 1440.97 200.00 19847.74 19647.74 2.66 7.15 509.00
14 4 HBACK 1 1264 2048.37 3291.24 783.58 1313.14 865.22 200.00 47296.61 47096.61 4.52 37.81 92.57
15 5 HDTOP 1 130 2267.78 5015.92 584.16 1064.98 484.53 200.00 32814.80 32614.80 4.26 19.88 439.93
16 6 MCARA 1 14 762.42 812.92 379.97 591.53 38.81 345.00 3230.60 2885.60 2.02 3.21 217.26
17 7 MIBUS 1 43 2700.11 4529.85 1286.59 1471.73 1068.48 200.00 20545.10 20345.10 2.67 6.29 690.80
18 8 PANVN 1 62 2146.99 3552.48 714.79 1333.94 626.92 200.00 22216.09 22016.09 3.41 14.50 451.17
19 9 RDSTR 1 2 684.73 685.51 684.73 684.73 718.66 200.00 1169.46 969.46 0.00 -2.75 484.73
110 10 SEDAN 1 1476 1816.82 2928.87 759.81 1144.24 829.98 200.00 29634.63 29434.63 4.09 24.18 76.24
111 11 STNWG 1 1173 2014.57 4063.93 734.19 1164.22 674.26 200.00 55922.13 55722.13 6.27 55.03 118.66
112 12 TRUCK 1 120 2662.47 4675.85 807.55 1392.42 672.78 200.00 22405.44 22205.44 2.81 7.30 426.84
113 13 UTE 1 260 2296.96 3728.77 782.51 1384.33 755.53 200.00 28012.83 27812.83 3.35 14.33 231.25
30
• Valores dos sinistros vs idade do motorista:
> with(dados2,describeBy(Quantia, Idademot, mat = TRUE,digits=2))
item group1 vars n mean sd median trimmed mad min max range skew kurtosis se
11 1 1 1 496 2635.83 4320.54 994.66 1705.46 1178.16 200 46868.18 46668.18 4.37 29.93 194.00
12 2 2 1 932 2129.66 4106.12 760.19 1241.36 830.55 200 55922.13 55722.13 5.64 47.35 134.50
13 3 3 1 1113 1915.64 3065.13 743.77 1201.50 806.19 200 31974.77 31774.77 3.74 19.19 91.88
14 4 4 1 1104 1943.21 3503.93 750.47 1165.36 816.13 200 47296.61 47096.61 5.36 43.88 105.46
15 5 5 1 614 1728.68 2798.37 702.76 1054.16 745.40 200 22216.09 22016.09 3.68 16.95 112.93
16 6 6 1 365 1872.79 3405.27 852.32 1110.64 826.48 200 31243.67 31043.67 4.69 27.66 178.24
31
BUS
CONVT
COUPE
HBACK
HDTOP
MCARA
MIBUS
PANVN
RDSTR
SEDAN
STNWG
TRUCK
UTE
0
10
20
30
40
50
Valor dos sinistros x $10.000
1
2
3
4
0
10
20
30
40
50
Idade d
o m
oto
rista
Valor dos sinistrosx $10.000
F
M
0
10
20
30
40
50
Sexo
Valor dos sinistrosx $10.000
A
B
C
D
E
F
0
10
20
30
40
50
Áre
a d
e re
sid
ência
Valor dos sinistrosx $10.000
1
2
3
4
5
6
0
10
20
30
40
50
Idade d
o m
oto
rista
Valor dos sinistrosx $10.000
0
2
4
6
8
10
12
14
0
10
20
30
40
50
Valo
r do ve
ículo
Valor dos sinistrosx $10.000
Fig
ur
a 7
– D
istribu
ição
do
s va
lores d
os sin
istros seg
un
do
as co
va
riáv
eis (tod
as a
s ap
ólices co
m sin
istro).
32
Exemplo 6 – Análise do desempenho de cinco tipos de turbinas de aviões. Foi conduzido um
experimento com 10 turbinas de cada tipo
Variável explicativa – Tipo de turbina, com níveis I, II, III, IV e V.
Variável resposta – Tempo de vida da turbina, em milhões de ciclos até verificada a perda de
velocidade.
Objetivos –
• Estimar parâmetros correspondentes às distribuições dos tempos de vida dos cinco tipos de
turbinas;
• Comparar os tempos médios de vida, identificar quais turbinas são mais resistentes.
33
Dados:
> turbdata
I II III IV V
1 3.03 3.19 3.46 5.88 6.43
2 5.53 4.26 5.22 6.74 9.97
3 5.60 4.47 5.69 6.90 10.39
4 9.30 4.53 6.54 6.98 13.55
5 9.92 4.67 9.16 7.21 14.45
6 12.51 4.69 9.40 8.14 14.72
7 12.95 5.78 10.19 8.59 16.81
8 15.21 6.79 10.71 9.80 18.39
9 16.04 9.37 12.58 12.28 20.84
10 16.84 12.75 13.41 25.46 21.51
34
Análise descritiva:
Turbina 1 Turbina 2 Turbina 3 Turbina 4 Turbina 5
510
1520
25
Tipo de turbina
Tem
po a
té p
erda
de
velo
cida
de (m
ilhõe
s de
cic
los)
Figura – Distribuição dos tempos de vida segundo tipo de turbina.
35
> medias=with(turbdata,tapply(Tempo,Turbina,mean)); medias
Turbina 1 Turbina 2 Turbina 3 Turbina 4 Turbina 5
10.693 6.050 8.636 9.798 14.706
> variancias=with(turbdata,tapply(Tempo,Turbina,var )); variancias
Turbina 1 Turbina 2 Turbina 3 Turbina 4 Turbina 5
23.225512 8.497489 10.828116 33.711796 23.652316
> cvs=sqrt(variancias)/medias;cvs # Coeficientes de variação.
Turbina 1 Turbina 2 Turbina 3 Turbina 4 Turbina 5
0.4506954 0.4818257 0.3810341 0.5925889 0.3307061