# Pacotes
library(plyr)
library(plotrix)Aplicação
A seguinte tabela foi reproduzida de PAULA, G.A. (2013) (https://www.ime.usp.br/~giapaula/textoregressao.htm), a qual resume os resultados de um estudo de seguimento em que médicos britânicos foram acompanhados durante a década de 50 e observada, em particular, a ocorrência de mortes por câncer de pulmão segundo o consumo médio diário de cigarros e a faixa etária.
Nossa análise se propõe a:
- Apresentar um modelo para descrever a taxa média de mortes
- Fazer um diagnóstico do modelo ajustado
- Apresentar uma estimativa para o melhor cenário
## Envelope e diagnósticos para Poisson baixados do site de Frederico Poleto
## http://www.poleto.com/funcoes.html
## Scripts sugeridos pelo Professor Caio Azevedo (IME- Unicamp)
## Exercício proposto pelo Professor Vanderly Janeiro (BPE - UEM)
source("envel_pois.R")
source("diag_pois.R")
dados <- read.table("breslow.dat")
tdados <- rbind(dados[1:4,1],dados[1:4,2],
dados[5:8,1],dados[5:8,2],
dados[9:12,1],dados[9:12,2],
dados[13:16,1],dados[13:16,2])
tdados <- matrix(tdados,nrow=8,byrow=FALSE)
dimnames(tdados)<- list(c("0","","1-9","","10-30","","+30",""),
c("40-49","50-59","60-69","70-80"))
names(dimnames(tdados)) <- c("Consumo médio diário de cigarros"," Faixa etária")
tdados## Faixa etária
## Consumo médio diário de cigarros 40-49 50-59 60-69 70-80
## 0 0.0 3.0 0.0 3.0
## 33679.0 21131.5 10599.0 4495.5
## 1-9 0.0 1.0 3.0 3.0
## 6002.5 4396.0 2813.5 1664.5
## 10-30 7.0 29.0 41.0 45.0
## 34414.5 25429.0 13271.0 4765.5
## +30 3.0 16.0 36.0 11.0
## 5881.0 6493.5 3466.5 769.0a. Análise descritiva.
Vamos estratificar o consumo de cigarros e faixas etárias como variáveis categóricas.
auxncig <- as.factor(dados[,3])
auxfe <- as.factor(dados[,4])
auxncig<-revalue(auxncig, c("1"="0", "2"="1-9","3"="10-30","4"="+30"))
auxfe<-revalue(auxfe, c("1"="40-49", "2"="50-59","3"="60-69","4"="70-80"))
dadosc <- data.frame(dados[,1],dados[,2],auxncig,auxfe)
colnames(dadosc) <- c("nmortes","panos","consumo","fetaria")
dadosc <- as.data.frame(dadosc)
dadosc## nmortes panos consumo fetaria
## 1 0 33679.0 0 40-49
## 2 3 21131.5 0 50-59
## 3 0 10599.0 0 60-69
## 4 3 4495.5 0 70-80
## 5 0 6002.5 1-9 40-49
## 6 1 4396.0 1-9 50-59
## 7 3 2813.5 1-9 60-69
## 8 3 1664.5 1-9 70-80
## 9 7 34414.5 10-30 40-49
## 10 29 25429.0 10-30 50-59
## 11 41 13271.0 10-30 60-69
## 12 45 4765.5 10-30 70-80
## 13 3 5881.0 +30 40-49
## 14 16 6493.5 +30 50-59
## 15 36 3466.5 +30 60-69
## 16 11 769.0 +30 70-80Temos uma array (matriz tridimensional) 4x4x4. Agora, vamos observar as taxas de mortes observadas
# Gráficos das taxas observadas
taxas <- dadosc$nmortes/dadosc$panos
par(mfrow=c(1,1))
plot(taxas[dadosc$consumo=="0"],pch=19,lty=1,type="b",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,ylim=c(0,0.01430429),xlab="faixa etária",ylab="taxa de mortes")
lines(taxas[dadosc$consumo=="1-9"],pch=17,lty=2,type="b",cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(taxas[dadosc$consumo=="10-30"],pch=15,lty=3,type="b",cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(taxas[dadosc$consumo=="+30"],pch=18,lty=4,type="b",cex=1.2,cex.axis=1.2,cex.lab=1.2)
axis(2,seq(0,0.015,length.out=7))
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
legend(1,0.015,c("0","1-9","10-30","+30"),lty=c(1,2,3,4),pch=c(19,17,15,18),bty="n",cex=1.3)
Vamos analisar o ajuste do modelo saturado
result<-fit.model <- glm(nmortes~offset(log(panos))+consumo+fetaria+consumo*fetaria,family=poisson(link="log"),data=dadosc)
summary(result)##
## Call:
## glm(formula = nmortes ~ offset(log(panos)) + consumo + fetaria +
## consumo * fetaria, family = poisson(link = "log"), data = dadosc)
##
## Deviance Residuals:
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.373e+01 6.965e+04 0 1
## consumo1-9 1.725e+00 9.851e+04 0 1
## consumo10-30 2.523e+01 6.965e+04 0 1
## consumo+30 2.615e+01 6.965e+04 0 1
## fetaria50-59 2.487e+01 6.965e+04 0 1
## fetaria60-69 1.156e+00 9.851e+04 0 1
## fetaria70-80 2.641e+01 6.965e+04 0 1
## consumo1-9:fetaria50-59 -1.253e+00 9.851e+04 0 1
## consumo10-30:fetaria50-59 -2.314e+01 6.965e+04 0 1
## consumo+30:fetaria50-59 -2.329e+01 6.965e+04 0 1
## consumo1-9:fetaria60-69 2.400e+01 1.206e+05 0 1
## consumo10-30:fetaria60-69 1.564e+00 9.851e+04 0 1
## consumo+30:fetaria60-69 1.857e+00 9.851e+04 0 1
## consumo1-9:fetaria70-80 -7.311e-01 9.851e+04 0 1
## consumo10-30:fetaria70-80 -2.258e+01 6.965e+04 0 1
## consumo+30:fetaria70-80 -2.308e+01 6.965e+04 0 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4.7258e+02 on 15 degrees of freedom
## Residual deviance: 4.5496e-10 on 0 degrees of freedom
## AIC: 83.479
##
## Number of Fisher Scoring iterations: 21 result##
## Call: glm(formula = nmortes ~ offset(log(panos)) + consumo + fetaria +
## consumo * fetaria, family = poisson(link = "log"), data = dadosc)
##
## Coefficients:
## (Intercept) consumo1-9
## -33.7272 1.7247
## consumo10-30 consumo+30
## 25.2269 26.1463
## fetaria50-59 fetaria60-69
## 24.8673 1.1561
## fetaria70-80 consumo1-9:fetaria50-59
## 26.4150 -1.2532
## consumo10-30:fetaria50-59 consumo+30:fetaria50-59
## -23.1433 -23.2924
## consumo1-9:fetaria60-69 consumo10-30:fetaria60-69
## 24.0028 1.5644
## consumo+30:fetaria60-69 consumo1-9:fetaria70-80
## 1.8574 -0.7311
## consumo10-30:fetaria70-80 consumo+30:fetaria70-80
## -22.5772 -23.0813
##
## Degrees of Freedom: 15 Total (i.e. Null); 0 Residual
## Null Deviance: 472.6
## Residual Deviance: 4.55e-10 AIC: 83.48 AICM0 <- AIC(fit.model)
BICM0 <- BIC(fit.model)
ez<-qnorm(0.975)
rebeta <- (summary(result))$coef
rebetaM1 <- cbind(rebeta[,1],rebeta[,2],rebeta[,1]-ez*rebeta[,2],rebeta[,1]+ez*rebeta[,2],rebeta[,3],rebeta[,4])
rownames(rebetaM1)<-c("alpha","beta_2","beta_3","beta_4",
"gamma_2","gamma_3","gamma_4",
"(alphabeta)_{22}","(alphabeta)_{23}",
"(alphabeta)_{24}","(alphabeta)_{32}",
"(alphabeta)_{33}","(alphabeta)_{34}",
"(alphabeta)_{42}","(alphabeta)_{43}",
"(alphabeta)_{44}")
rebetaM1## [,1] [,2] [,3] [,4] [,5]
## alpha -33.7272149 69653.80 -136552.7 136485.2 -4.842121e-04
## beta_2 1.7246984 98505.35 -193065.2 193068.7 1.750868e-05
## beta_3 25.2268918 69653.80 -136493.7 136544.2 3.621754e-04
## beta_4 26.1463451 69653.80 -136492.8 136545.1 3.753757e-04
## gamma_2 24.8673071 69653.80 -136494.1 136543.8 3.570129e-04
## gamma_3 1.1561148 98505.35 -193065.8 193068.1 1.173657e-05
## gamma_4 26.4149950 69653.80 -136492.5 136545.4 3.792326e-04
## (alphabeta)_{22} -1.2532410 98505.35 -193068.2 193065.7 -1.272257e-05
## (alphabeta)_{23} -23.1433336 69653.80 -136542.1 136495.8 -3.322623e-04
## (alphabeta)_{24} -23.2924055 69653.80 -136542.2 136495.6 -3.344025e-04
## (alphabeta)_{32} 24.0028293 120643.92 -236433.7 236481.7 1.989560e-04
## (alphabeta)_{33} 1.5644439 98505.35 -193065.4 193068.5 1.588182e-05
## (alphabeta)_{34} 1.8573732 98505.35 -193065.1 193068.8 1.885556e-05
## (alphabeta)_{42} -0.7311463 98505.35 -193067.7 193066.2 -7.422402e-06
## (alphabeta)_{43} -22.5771671 69653.80 -136541.5 136496.4 -3.241340e-04
## (alphabeta)_{44} -23.0813209 69653.80 -136542.0 136495.9 -3.313720e-04
## [,6]
## alpha 0.9996137
## beta_2 0.9999860
## beta_3 0.9997110
## beta_4 0.9997005
## gamma_2 0.9997151
## gamma_3 0.9999906
## gamma_4 0.9996974
## (alphabeta)_{22} 0.9999898
## (alphabeta)_{23} 0.9997349
## (alphabeta)_{24} 0.9997332
## (alphabeta)_{32} 0.9998413
## (alphabeta)_{33} 0.9999873
## (alphabeta)_{34} 0.9999850
## (alphabeta)_{42} 0.9999941
## (alphabeta)_{43} 0.9997414
## (alphabeta)_{44} 0.9997356O modelo saturado alcançou AIC de 83.48, indicando que há fortes indícios de interação entre a faixa de consumo de cigarros e a faixa etária para explicar o número de mortes por câncer de pulmão. A seguir, vamos testar o modelo reduzido, sem considerar as interações, para verificar a perda de qualidade de ajuste e eventual Paradoxo de Simpson (confundidores).
b.
Vamos testar o modelo sem interação
# Modelo sem interação
result<-fit.model <- glm(nmortes~offset(log(panos))+consumo+fetaria,family=poisson(link="log"),data=dadosc)
summary(fit.model)##
## Call:
## glm(formula = nmortes ~ offset(log(panos)) + consumo + fetaria,
## family = poisson(link = "log"), data = dadosc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.04914 -0.74407 -0.03998 0.47262 1.21622
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.4241 0.5094 -22.424 < 2e-16 ***
## consumo1-9 1.4088 0.5566 2.531 0.0114 *
## consumo10-30 2.8662 0.4182 6.854 7.20e-12 ***
## consumo+30 3.7580 0.4269 8.803 < 2e-16 ***
## fetaria50-59 1.7693 0.3473 5.095 3.49e-07 ***
## fetaria60-69 2.8973 0.3357 8.630 < 2e-16 ***
## fetaria70-80 3.7914 0.3409 11.121 < 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: 472.580 on 15 degrees of freedom
## Residual deviance: 11.909 on 9 degrees of freedom
## AIC: 77.388
##
## Number of Fisher Scoring iterations: 5 summary(result)##
## Call:
## glm(formula = nmortes ~ offset(log(panos)) + consumo + fetaria,
## family = poisson(link = "log"), data = dadosc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.04914 -0.74407 -0.03998 0.47262 1.21622
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.4241 0.5094 -22.424 < 2e-16 ***
## consumo1-9 1.4088 0.5566 2.531 0.0114 *
## consumo10-30 2.8662 0.4182 6.854 7.20e-12 ***
## consumo+30 3.7580 0.4269 8.803 < 2e-16 ***
## fetaria50-59 1.7693 0.3473 5.095 3.49e-07 ***
## fetaria60-69 2.8973 0.3357 8.630 < 2e-16 ***
## fetaria70-80 3.7914 0.3409 11.121 < 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: 472.580 on 15 degrees of freedom
## Residual deviance: 11.909 on 9 degrees of freedom
## AIC: 77.388
##
## Number of Fisher Scoring iterations: 5Este parece ser um bom modelo para o conjunto de dados, uma vez que o AIC é de 77.39 e todos os coeficientes possuem p-valor menor que 0.05. Tendo em vista que é mais simples que o modelo saturado e a perda de qualidade do ajuste é aceitável, apontamos como sendo o melhor modelo:
\(y=-11.42+1.41x_{1}+2.87x_{2}+3.76x_{3}+1.77x_{4}+2.90x_{5}+3.79x_{6}\)
Em que
\(y>0\) é o número de mortes esperadas por câncer de pulmão;
\(x_{1}=1\) caso o consumo seja de 1 a 9 cigarros e \(x_{1}=0\) caso contrário;
\(x_{2}=1\) caso o consumo seja de 10 a 30 cigarros e \(x_{2}=0\) caso contrário;
\(x_{3}=1\) caso o consumo seja maior que 30 cigarros e \(x_{3}=0\) caso contrário;
\(x_{4}=1\) caso a idade seja entre 50 e 59 anos e \(x_{4}=0\) caso contrário;
\(x_{5}=1\) caso a idade seja entre 60 e 69 anos e \(x_{5}=0\) caso contrário;
\(x_{6}=1\) caso a idade seja entre 70 e 80 anos e \(x_{6}=0\) caso contrário;
b. Diagnóstico.
Vamos ver a simulação de envelope e realizar o teste diagnóstico para o Modelo Log-Poisson:
envel.pois(fit.model)
## Banda de 90 % de confianca, obtida por 100 simulacoes.diag.pois(fit.model)
## $ResQuantil
## [1] 0.293387629 1.186526455 -1.198265241 0.331540265 0.334657010
## [6] -0.081893733 0.457411758 -0.236004751 0.066802937 0.006691269
## [11] -0.792983735 0.664021396 0.260780566 -0.368294654 1.253022496
## [16] -1.178531617
##
## $ResCompDesv
## 1 2 3 4 5 6
## -0.81599007 1.05473157 -1.62831874 0.41354537 -0.70857250 -0.13155983
## 7 8 9 10 11 12
## 0.36685029 -0.11885468 0.08573159 0.03676170 -0.38799866 0.33595728
## 13 14 15 16
## 0.12095993 -0.30524049 0.70499176 -0.96893244
##
## $ResAnscombe
## 1 2 3 4 5 6
## -0.91004892 1.22014282 -2.17344568 0.52710931 -0.77709122 -0.14592804
## 7 8 9 10 11 12
## 0.45470630 -0.16693871 0.15083465 0.06595474 -0.77838164 0.68541692
## 13 14 15 16
## 0.14547241 -0.44482744 1.16944374 -1.31829960
##
## $ResPearsonStd
## 1 2 3 4 5 6 7
## -0.6379360 1.6297554 -1.8234365 0.7100452 -0.5356637 -0.1581274 0.5907944
## 8 9 10 11 12 13 14
## -0.2308366 0.2679502 0.1185730 -1.5312489 1.4231007 0.1774775 -0.6366997
## 15 16
## 2.0075993 -1.6900282
##
## $ResWilliams
## 1 2 3 4 5 6 7
## -0.8006925 1.2227708 -1.7028344 0.5467199 -0.6986899 -0.1369226 0.4575428
## 8 9 10 11 12 13 14
## -0.1828505 0.2257797 0.1005561 -1.3414407 1.2512923 0.1408398 -0.5082892
## 15 16
## 1.6570560 -1.3487623
##
## $Di
## 1 2 3 4 5 6
## 0.0061406764 0.1250859730 0.2772392713 0.0448723277 0.0028329255 0.0008220323
## 7 8 9 10 11 12
## 0.0266874206 0.0074038767 0.0214913008 0.0044565734 1.0129532197 0.9148060586
## 13 14 15 16
## 0.0020080943 0.0650643578 1.0079434715 0.3463709847
##
## $Ci
## 1 2 3 4 5 6 7
## 0.30070257 0.68666525 1.41057899 0.37012438 0.21121800 0.07156189 0.30431757
## 8 9 10 11 12 13 14
## 0.13291149 0.14071469 0.06209139 0.76506862 0.67738285 0.09162448 0.36685995
## 15 16
## 1.05766066 1.01225881
##
## $Dmax
## [1] 0.1049169284 0.4897588807 0.6376799675 0.2247860883 0.0004963952
## [6] 0.0016560143 0.0121281512 0.0040509076 0.0024606629 0.0022572937
## [11] 0.0643779581 0.0799926047 0.0031895335 0.0006545351 0.0372603227
## [16] 0.1228809863
##
## $h
## [1] 0.09553288 0.24792612 0.36855754 0.38386607 0.06464365 0.18707753
## [7] 0.34862731 0.49306277 0.67693279 0.68932969 0.75149710 0.75972836
## [13] 0.30856511 0.52907865 0.63643935 0.45913508Medidas de AIC e BIC do modelo:
AICM1 <- AIC(fit.model)
AICM1## [1] 77.38829BICM1 <- BIC(fit.model)
BICM1## [1] 82.79641ez<-qnorm(0.975)
rebeta <- (summary(result))$coef
rebeta## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.424074 0.5094497 -22.424343 2.278338e-111
## consumo1-9 1.408817 0.5566427 2.530918 1.137644e-02
## consumo10-30 2.866195 0.4182040 6.853580 7.202441e-12
## consumo+30 3.758026 0.4269115 8.802821 1.334187e-18
## fetaria50-59 1.769257 0.3472645 5.094839 3.490369e-07
## fetaria60-69 2.897256 0.3357195 8.629992 6.135653e-18
## fetaria70-80 3.791443 0.3409321 11.120817 9.935772e-29Matriz de Covariância:
rebetaM1 <- cbind(rebeta[,1],rebeta[,2],rebeta[,1]-ez*rebeta[,2],rebeta[,1]+ez*rebeta[,2],rebeta[,3],rebeta[,4])
mrebeta <- rbind(rebetaM1)
covbetaM1 <- vcov(fit.model)
covbetaM1## (Intercept) consumo1-9 consumo10-30 consumo+30 fetaria50-59
## (Intercept) 0.25953895 -0.164265127 -0.1659998098 -1.648962e-01 -0.0985718718
## consumo1-9 -0.16426513 0.309851149 0.1666076526 1.664102e-01 -0.0010347237
## consumo10-30 -0.16599981 0.166607653 0.1748945953 1.667832e-01 -0.0009096633
## consumo+30 -0.16489623 0.166410245 0.1667831647 1.822535e-01 -0.0029007861
## fetaria50-59 -0.09857187 -0.001034724 -0.0009096633 -2.900786e-03 0.1205926497
## fetaria60-69 -0.09835108 -0.001979938 -0.0010844152 -3.190943e-03 0.1001937798
## fetaria70-80 -0.09974588 -0.004060856 -0.0002223685 6.723625e-06 0.0999689983
## fetaria60-69 fetaria70-80
## (Intercept) -0.098351079 -9.974588e-02
## consumo1-9 -0.001979938 -4.060856e-03
## consumo10-30 -0.001084415 -2.223685e-04
## consumo+30 -0.003190943 6.723625e-06
## fetaria50-59 0.100193780 9.996900e-02
## fetaria60-69 0.112707568 9.998753e-02
## fetaria70-80 0.099987529 1.162347e-01c. Estimativa.
Gráfico com as estimativas para cada \(\beta\) e \(\delta\):
par(mfrow=c(1,2))
plotCI(rebetaM1[2:4,1],li=rebetaM1[2:4,3],ui=rebetaM1[2:4,4],pch=19,cex=1.2,cex.lab=1.2,cex.main=1.2,xlab="parâmetro",ylab="estimativa",axes=F)
axis(2,cex=1.2,at=seq(0,5,by=1),cex.lab=1.2,cex.main=1.2,cex.axis=1.2)
axis(1,c(1,2,3),labels=c(expression(beta[2]),expression(beta[3]),expression(beta[4])),cex=1.2,cex.lab=1.2,cex.main=1.2,cex.axis=1.2)
plotCI(rebetaM1[5:7,1],li=rebetaM1[5:7,3],ui=rebetaM1[5:7,4],pch=19,cex=1.2,cex.lab=1.2,cex.main=1.2,xlab="parâmetro",ylab="estimativa",axes=F)
axis(2,cex=1.2,at=seq(0,5,by=1),cex.lab=1.2,cex.main=1.2)
axis(1,c(1,2,3),labels=c(expression(delta[2]),expression(delta[3]),expression(delta[4])),cex=1.2,cex.lab=1.2,cex.main=1.2,cex.axis=1.2)
Vamos comparar os dados observados com o modelo preditivo.
pred<-predict(fit.model,type=c("response"),se.fit = TRUE)
mupredM1 <- unique(pred$fit)
semupredM1 <- unique(pred$se.fit)
liICM1 = mupredM1-ez*semupredM1
lsICM1 = mupredM1+ez*semupredM1
par(mfrow=c(2,2))
plot(dadosc$nmortes[dadosc$consumo=="0"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,4),main="consumo 0")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(mupredM1[dadosc$consumo=="0"],li=liICM1[dadosc$consumo=="0"],ui=lsICM1[dadosc$consumo=="0"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(mupredM1[dadosc$consumo=="0"],lwd=2)
plot(dadosc$nmortes[dadosc$consumo=="1-9"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,7),main="consumo 1-9")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(mupredM1[dadosc$consumo=="1-9"],li=liICM1[dadosc$consumo=="1-9"],ui=lsICM1[dadosc$consumo=="1-9"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(mupredM1[dadosc$consumo=="1-9"],lwd=2)
plot(dadosc$nmortes[dadosc$consumo=="10-30"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(2,60),main="consumo 10-30")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(mupredM1[dadosc$consumo=="10-30"],li=liICM1[dadosc$consumo=="10-30"],ui=lsICM1[dadosc$consumo=="10-30"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(mupredM1[dadosc$consumo=="10-30"],lwd=2)
plot(dadosc$nmortes[dadosc$consumo=="+30"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,45),main="consumo +30")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(mupredM1[dadosc$consumo=="+30"],li=liICM1[dadosc$consumo=="+30"],ui=lsICM1[dadosc$consumo=="+30"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(mupredM1[dadosc$consumo=="+30"],lwd=2)
Taxas de mortes preditas:
# Taxas de morte preditas
taxapredM1 <- mupredM1/dadosc$panos
setaxapredM1 <- semupredM1/dadosc$panos
liICM1 = taxapredM1-ez*setaxapredM1
lsICM1 = taxapredM1+ez*setaxapredM1
par(mfrow=c(2,2))
plot(taxas[dadosc$consumo=="0"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,0.001),main="consumo 0")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(taxapredM1[dadosc$consumo=="0"],li=liICM1[dadosc$consumo=="0"],ui=lsICM1[dadosc$consumo=="0"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(taxapredM1[dadosc$consumo=="0"],lwd=2)
plot(taxas[dadosc$consumo=="1-9"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,0.004),main="consumo 1-9")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(taxapredM1[dadosc$consumo=="1-9"],li=liICM1[dadosc$consumo=="1-9"],ui=lsICM1[dadosc$consumo=="1-9"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(taxapredM1[dadosc$consumo=="1-9"],lwd=2)
plot(taxas[dadosc$consumo=="10-30"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,0.015),main="consumo 10-30")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(taxapredM1[dadosc$consumo=="10-30"],li=liICM1[dadosc$consumo=="10-30"],ui=lsICM1[dadosc$consumo=="10-30"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(taxapredM1[dadosc$consumo=="10-30"],lwd=2)
plot(taxas[dadosc$consumo=="+30"],pch=1,lty=1,type="p",cex=1.2,cex.axis=1.2,cex.lab=1.2,axes=F,xlab="faixa etária",ylab="taxa de mortes",ylim=c(0,0.03),main="consumo +30")
axis(2)
axis(1,c(1,2,3,4),labels=c("40-49","50-59","60-69","70-80"))
plotCI(taxapredM1[dadosc$consumo=="+30"],li=liICM1[dadosc$consumo=="+30"],ui=lsICM1[dadosc$consumo=="+30"],add=TRUE,pch=19,cex=1.2,cex.axis=1.2,cex.lab=1.2)
lines(taxapredM1[dadosc$consumo=="+30"],lwd=2)
Conclusão:
Mesmo havendo uma considerável discrepância entre o observado e o predito nas faixas etárias mais altas, o modelo proposto é subsistente. Devemos ter em mente que o modelo preditivo não leva em consideração a morte em decorrência de outras causas concorrentes com o tabagismo. Em especial, na faixa etária entre 70 e 80 anos estamos além do limiar da expectativa de vida ao nascer, portanto, é esperado um decréscimo na mortalidade pelo tabagismo nessa população.
Compare nossa análise e resultados com a página 309 da obra supracitada.
* www.linkedin.com/in/msrcos3s