Dúvida no ajuste de modelo de GEE para dados binários, longitudinais
Prezados(as) Ajustei o seguinte modelo de GEE para dados binários, longitudinais: load(https://dl.dropboxusercontent.com/u/73337918/d6.rda) attach(d6) d6$REF_P <- as.factor(d6$REF_P) d6$TCM_PRE <- as.factor(d6$TCM_PRE) d6$TCM_REL <- as.factor(d6$TCM_REL) d6$MAN_CON <- as.factor(d6$MAN_CON) d6$MAN_PER <- as.factor(d6$MAN_PER) d6$EST_CIV <- as.factor(d6$EST_CIV) d6$GR_ESC <- as.factor(d6$GR_ESC) d6$GR_OCUP <- as.factor(d6$GR_OCUP) d6$GR_ESCP <- as.factor(d6$GR_ESCP) d6$INT_PGOV <- as.factor(d6$INT_PGOV) d6$INT_PPRE <- as.factor(d6$INT_PPRE) d6$INT_PART <- as.factor(d6$INT_PART) d6$GR_DCAP <- as.factor(d6$GR_DCAP) d6$RESC_ANT <- as.factor(d6$RESC_ANT) # d6$GR_ESC <- relevel(d6$GR_ESC, ref="LE_FI") d6$REG_MES <- relevel(d6$REG_MES, ref="Ms_Mt") library(geepack) mod6c <- geeglm(RES_CONT ~ GR_ESC+GR_OCUP+SEXO+MAN_PER+REF_P+ TCM_PRE+TCM_REL+INT_PGOV+INT_PPRE+REG_MES+GR_DCAP+IDHM_PNUD+RESC_ANT,id = MUNIC, data = d6, family = binomial(link=logit), corstr = "ar1") summary(mod6c) Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) 1.54944 0.58054 7.123 0.00761 ** GR_ESCFC_MI 0.05587 0.15362 0.132 0.71610 GR_ESCMC_SI 0.25302 0.12323 4.215 0.04006 * GR_ESCSC_PG 0.47851 0.15347 9.722 0.00182 ** GR_OCUP2 -0.01318 0.19240 0.005 0.94540 GR_OCUP3 -0.49723 0.20344 5.974 0.01452 * GR_OCUP4 0.20693 0.18915 1.197 0.27397 GR_OCUP5 0.05465 0.13705 0.159 0.69005 GR_OCUP6 -0.16555 0.20084 0.679 0.40979 GR_OCUP7 0.25329 0.22578 1.259 0.26192 GR_OCUP8 -0.56514 0.20535 7.574 0.00592 ** GR_OCUP9 0.11389 0.23092 0.243 0.62187 GR_OCUP10 -0.32457 0.16519 3.860 0.04944 * SEXOM -0.38449 0.15291 6.322 0.01192 * MAN_PER2 -0.02760 0.08483 0.106 0.74493 MAN_PER3 0.89274 0.22867 15.241 9.46e-05 *** REF_P2 -0.90280 0.16081 31.516 1.98e-08 *** TCM_PRE2 0.32531 0.17539 3.440 0.06362 . TCM_PRE3 1.10453 0.38724 8.136 0.00434 ** TCM_PRE4 0.07393 0.16962 0.190 0.66292 TCM_PRE5 0.45656 0.17748 6.618 0.01010 * TCM_PRE6 -1.07012 0.43169 6.145 0.01318 * TCM_REL2 -0.88695 0.22036 16.202 5.69e-05 *** TCM_REL3 -0.51005 0.23513 4.706 0.03006 * TCM_REL4 -0.20505 0.21699 0.893 0.34466 TCM_REL5 0.05124 0.26072 0.039 0.84418 TCM_REL6 -0.53652 0.27731 3.743 0.05303 . TCM_REL7 -0.49148 0.21464 5.243 0.02204 * TCM_REL8 0.21678 0.21973 0.973 0.32385 TCM_REL9 0.32581 0.31959 1.039 0.30799 TCM_REL10 0.59204 0.25881 5.233 0.02216 * TCM_REL11 -0.06144 0.22802 0.073 0.78759 TCM_REL12 -0.51748 0.28703 3.250 0.07140 . INT_PGOV1 -0.86442 0.37774 5.237 0.02211 * INT_PGOV2 -0.62852 0.39466 2.536 0.11126 INT_PPRE1 1.07552 0.37019 8.441 0.00367 ** INT_PPRE2 0.94990 0.43197 4.836 0.02788 * REG_MESMs_CN -0.29830 0.19545 2.329 0.12696 REG_MESMs_CS -0.25052 0.22930 1.194 0.27460 REG_MESMs_EO 0.01757 0.31662 0.003 0.95575 REG_MESMs_Nd -0.41248 0.20274 4.139 0.04190 * REG_MESMs_S -1.09423 0.22003 24.732 6.59e-07 *** REG_MESMs_VSF -0.68058 0.26449 6.621 0.01008 * GR_DCAP2 0.05559 0.15406 0.130 0.71821 GR_DCAP3 -0.02412 0.18018 0.018 0.89351 GR_DCAP4 0.21081 0.21921 0.925 0.33622 GR_DCAP5 0.62488 0.22953 7.412 0.00648 ** GR_DCAP6 0.60297 0.74007 0.664 0.41522 IDHM_PNUD -2.48469 0.77384 10.310 0.00132 ** RESC_ANT1 1.60297 0.08113 390.413 < 2e-16 *** Não consegui entender os resultados seguintes: INT_PGOV1 -0.86442 0.37774 5.237 0.02211 * INT_PGOV2 -0.62852 0.39466 2.536 0.11126 ..............................................
tapply(fitted(mod6c), INT_PGOV, mean) 0 1 2 0.6868 0.7403 0.7717 tapply(predict(mod6c), INT_PGOV, mean) 0 1 2 0.9619 1.2828 1.4899
A variável categórica INT_PGOV tem três níveis: 0, 1, e 2. O R coloca o 0 como baseline. Se o preditor linear e o valor ajustado da resposta são maiores nos níveis 1 e 2 de INT_PGOV, as estimativas dos parâmetros correspondentes a (INT_PGOV1:-0.86442) e (INT_PGOV2: -0.62852) não deviam ser positivas ????. Alguém tem alguma ideia? Desde já agradeço qualquer comentário. Att. -- Gilênio Borges Fernandes Professor Associado IV (Aposentado) Professor Adjunto A (Substituto) Universidade Federal da Bahia Instituto de Matemática Departamento de Estatística Av. Adhemar de Barros, s/n – Ondina. 40.170-110 - Salvador - BA, Brasil Tel.: (071)3283-6340/6341/6337 Fax: (071)3283-6336 URL: http://lattes.cnpq.br/6764860618464860
On 11/08/15 14:50, Gilenio Borges Fernandes wrote:
Se o preditor linear e o valor ajustado da resposta são maiores nos níveis 1 e 2 de INT_PGOV, as estimativas dos parâmetros correspondentes a (INT_PGOV1:-0.86442) e (INT_PGOV2: -0.62852) não deviam ser positivas ????.
com mais de uma covariavel essa afirmativa e' apenas valida na condicao de ortogonalidade.
Isso e' so' um dos problemas... ha' mais. Veja tres problemas nesse exemplo: ### possiveis efeitos de colinearidade ### covariaveis dat <- data.frame(a=gl(3, 50)) n <- nrow(dat) dat$xc <- rnorm(n, unclass(dat$a), 0.3) ### nao ortogonal dat$xo <- rnorm(n, 1:3, 0.3) ### ortogonal cor(model.matrix(~0+a, dat), dat[c('xc', 'xo')]) ### simula do modelo linear beta <- c(5, 1:2, -2) dat$yc <- rnorm(n, model.matrix(~a+xc,dat)%*%beta, 0.2) dat$yo <- rnorm(n, model.matrix(~a+xo,dat)%*%beta, 0.2) ### 1. condicionalidade versus marginalidade ### -> medias por grupo nao indicam efeito condicional sapply(dat[c('yc', 'yo')], tapply, dat['a'], mean) ### estimativas dos coeficientes mc <- lm(yc~a+xc, dat) mo <- lm(yo~a+xo, dat) ### 2. estimativas viciadas cbind(colin=coef(mc), ortog=coef(mo)) - beta ### 3. inflacao de variancia diag(vcov(mc))/diag(vcov(mo)) ### tamanho do prejuizo! Elias On 11/08/15 14:56, Elias Teixeira Krainski wrote:
On 11/08/15 14:50, Gilenio Borges Fernandes wrote:
Se o preditor linear e o valor ajustado da resposta são maiores nos níveis 1 e 2 de INT_PGOV, as estimativas dos parâmetros correspondentes a (INT_PGOV1:-0.86442) e (INT_PGOV2: -0.62852) não deviam ser positivas ????.
com mais de uma covariavel essa afirmativa e' apenas valida na condicao de ortogonalidade. _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Eu não consegui rodar este script, tem alguma biblioteca para isto? Prezados(as) Ajustei o seguinte modelo de GEE para dados binários, longitudinais: load(https://dl.dropboxusercontent.com/u/73337918/d6.rda) attach(d6) d6$REF_P <- as.factor(d6$REF_P) d6$TCM_PRE <- as.factor(d6$TCM_PRE) d6$TCM_REL <- as.factor(d6$TCM_REL) d6$MAN_CON <- as.factor(d6$MAN_CON) d6$MAN_PER <- as.factor(d6$MAN_PER) d6$EST_CIV <- as.factor(d6$EST_CIV) d6$GR_ESC <- as.factor(d6$GR_ESC) d6$GR_OCUP <- as.factor(d6$GR_OCUP) d6$GR_ESCP <- as.factor(d6$GR_ESCP) d6$INT_PGOV <- as.factor(d6$INT_PGOV) d6$INT_PPRE <- as.factor(d6$INT_PPRE) d6$INT_PART <- as.factor(d6$INT_PART) d6$GR_DCAP <- as.factor(d6$GR_DCAP) d6$RESC_ANT <- as.factor(d6$RESC_ANT) # d6$GR_ESC <- relevel(d6$GR_ESC, ref="LE_FI") d6$REG_MES <- relevel(d6$REG_MES, ref="Ms_Mt") library(geepack) mod6c <- geeglm(RES_CONT ~ GR_ESC+GR_OCUP+SEXO+MAN_PER+REF_P+ TCM_PRE+TCM_REL+INT_PGOV+INT_PPRE+REG_MES+GR_DCAP+IDHM_PNUD+RESC_ANT,id = MUNIC, data = d6, family = binomial(link=logit), corstr = "ar1") summary(mod6c) Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) 1.54944 0.58054 7.123 0.00761 ** GR_ESCFC_MI 0.05587 0.15362 0.132 0.71610 GR_ESCMC_SI 0.25302 0.12323 4.215 0.04006 * GR_ESCSC_PG 0.47851 0.15347 9.722 0.00182 ** GR_OCUP2 -0.01318 0.19240 0.005 0.94540 GR_OCUP3 -0.49723 0.20344 5.974 0.01452 * GR_OCUP4 0.20693 0.18915 1.197 0.27397 GR_OCUP5 0.05465 0.13705 0.159 0.69005 GR_OCUP6 -0.16555 0.20084 0.679 0.40979 GR_OCUP7 0.25329 0.22578 1.259 0.26192 GR_OCUP8 -0.56514 0.20535 7.574 0.00592 ** GR_OCUP9 0.11389 0.23092 0.243 0.62187 GR_OCUP10 -0.32457 0.16519 3.860 0.04944 * SEXOM -0.38449 0.15291 6.322 0.01192 * MAN_PER2 -0.02760 0.08483 0.106 0.74493 MAN_PER3 0.89274 0.22867 15.241 9.46e-05 *** REF_P2 -0.90280 0.16081 31.516 1.98e-08 *** TCM_PRE2 0.32531 0.17539 3.440 0.06362 . TCM_PRE3 1.10453 0.38724 8.136 0.00434 ** TCM_PRE4 0.07393 0.16962 0.190 0.66292 TCM_PRE5 0.45656 0.17748 6.618 0.01010 * TCM_PRE6 -1.07012 0.43169 6.145 0.01318 * TCM_REL2 -0.88695 0.22036 16.202 5.69e-05 *** TCM_REL3 -0.51005 0.23513 4.706 0.03006 * TCM_REL4 -0.20505 0.21699 0.893 0.34466 TCM_REL5 0.05124 0.26072 0.039 0.84418 TCM_REL6 -0.53652 0.27731 3.743 0.05303 . TCM_REL7 -0.49148 0.21464 5.243 0.02204 * TCM_REL8 0.21678 0.21973 0.973 0.32385 TCM_REL9 0.32581 0.31959 1.039 0.30799 TCM_REL10 0.59204 0.25881 5.233 0.02216 * TCM_REL11 -0.06144 0.22802 0.073 0.78759 TCM_REL12 -0.51748 0.28703 3.250 0.07140 . INT_PGOV1 -0.86442 0.37774 5.237 0.02211 * INT_PGOV2 -0.62852 0.39466 2.536 0.11126 INT_PPRE1 1.07552 0.37019 8.441 0.00367 ** INT_PPRE2 0.94990 0.43197 4.836 0.02788 * REG_MESMs_CN -0.29830 0.19545 2.329 0.12696 REG_MESMs_CS -0.25052 0.22930 1.194 0.27460 REG_MESMs_EO 0.01757 0.31662 0.003 0.95575 REG_MESMs_Nd -0.41248 0.20274 4.139 0.04190 * REG_MESMs_S -1.09423 0.22003 24.732 6.59e-07 *** REG_MESMs_VSF -0.68058 0.26449 6.621 0.01008 * GR_DCAP2 0.05559 0.15406 0.130 0.71821 GR_DCAP3 -0.02412 0.18018 0.018 0.89351 GR_DCAP4 0.21081 0.21921 0.925 0.33622 GR_DCAP5 0.62488 0.22953 7.412 0.00648 ** GR_DCAP6 0.60297 0.74007 0.664 0.41522 IDHM_PNUD -2.48469 0.77384 10.310 0.00132 ** RESC_ANT1 1.60297 0.08113 390.413 < 2e-16 *** Não consegui entender os resultados seguintes: INT_PGOV1 -0.86442 0.37774 5.237 0.02211 * INT_PGOV2 -0.62852 0.39466 2.536 0.11126 ..............................................
tapply(fitted(mod6c), INT_PGOV, mean)
0 1 2 0.6868 0.7403 0.7717
tapply(predict(mod6c), INT_PGOV, mean)
0 1 2 0.9619 1.2828 1.4899 A variável categórica INT_PGOV tem três níveis: 0, 1, e 2. O R coloca o 0 como baseline. Se o preditor linear e o valor ajustado da resposta são maiores nos níveis 1 e 2 de INT_PGOV, as estimativas dos parâmetros correspondentes a (INT_PGOV1:-0.86442) e (INT_PGOV2: -0.62852) não deviam ser positivas ????. Alguém tem alguma ideia? Desde já agradeço qualquer comentário. Att. -- Gilênio Borges Fernandes Professor Associado IV (Aposentado) Professor Adjunto A (Substituto) Universidade Federal da Bahia Instituto de Matemática Departamento de Estatística Av. Adhemar de Barros, s/n – Ondina. 40.170-110 - Salvador - BA, Brasil Tel.: (071)3283-6340/6341/6337 Fax: (071)3283-6336 URL: http://lattes.cnpq.br/6764860618464860 --- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus
participantes (3)
-
Elias Teixeira Krainski -
Gilenio Borges Fernandes -
Mauro Sznelwar