########################################################################################## ## Analise com os dados do IQVT ########################################################### ########################################################################################## ## #load("resultadosIQVT.RData") dados <- read.table("dadosBeta.txt", header=TRUE, sep="", dec= ".") ########################################################################################## ## Analise Descritiva #################################################################### ########################################################################################## require(lattice) xyplot(GLOBAL ~ RENDA | ESTADO, data=dados) ########################################################################################## ## Modelo Linear Simples ################################################################# ########################################################################################## mod.lm <- lm(GLOBAL ~ RENDA + ESTADO, data=dados) ########################################################################################## ## Modelos Transformados ################################################################# ########################################################################################## source("functions.R") source("functionsKw.R") source("genericRegressao.R") source("bandafunctions.R") require(latticeExtra) require(VGAM) require(rootSolve) require(numDeriv) ########################################################################################## ## Modelos transformados ################################################################# ########################################################################################## ## Logistica dados$y <- dados$GLOBAL dados$cov <- log(dados$RENDA) - mean(log(dados$RENDA)) tra.logit <- lm(logit(y) ~ cov + ESTADO, data=dados) ll.logit <- sum(log(1/dados$y + 1/(1-dados$y))) + logLik(tra.logit) ## Probit tra.probit <- lm(probit(y) ~ cov + ESTADO, data=dados) summary(tra.probit) ll.probit <- sum(log(diag(gradient(qnorm, x = as.numeric(dados$y))))) + logLik(tra.probit) ## Complemento Log-Log tra.cloglog <- lm(cloglog(y) ~ cov + ESTADO, data=dados) summary(tra.cloglog) ll.cloglog <- sum( log( - 1 / (log(1-dados$y)*(1-dados$y)))) + logLik(tra.cloglog) ## Loglog tra.loglog <- lm(loglog(y) ~ cov + ESTADO, data=dados) summary(tra.loglog) ll.loglog <- sum(log( - 1/(dados$y*log(dados$y)))) + logLik(tra.loglog) ## Cauchit tra.cauchit <- lm(cauchit(y) ~ cov + ESTADO, data=dados ) summary(tra.cauchit) argum <- pi*(dados$y - 0.5) ll.cauchit <- sum(log(pi*(1/cos(argum))^2)) + logLik(tra.cauchit) ## Este modelo está estranho alguma coisa estranha mod.ordaz <- lm(y ~ cov + ESTADO, data=dados) ## Escolhendo o lambda ll.aranda <- aranda.ordaz(mod.ordaz, lambda = c(0.01,10)) tra.ordaz <- lm(aranda1(y, lambda = 9.18) ~ RENDA + PORTE + ESTADO, data=dados) ll.aranda ########################################################################################## ## Base de dados para a construção das bandas de confiança ############################### ########################################################################################## ESTADO <- levels(dados$ESTADO) cov <- seq(-1,1.5, l= 100) newdados <- expand.grid(cov, ESTADO) names(newdados) <- c("cov", "ESTADO") ########################################################################################## ## Modelos Beta ########################################################################## ########################################################################################## require(bbmle) ## Link logit inits <- coef(tra.logit) inits <- coef(pp.beta.logit) beta.logit <- genericModel(formula = "y ~ cov + ESTADO", data = dados, inits = c(inits), family="beta", link = "logit") pp.beta.logit <- profile(beta.logit, tol.newmin=1) plot(pp.beta.logit) ## Plotando as bandas de confiança int.beta.logit <- cria.intervalo(beta.logit, formu= "~cov + ESTADO", data=newdados, nuisance=11, condicional=FALSE) ## Plotando o grafico p1 <- xyplot(y ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(inv.logit(preditos) + inv.logit(lim.inf) + inv.logit(lim.sup) ~ cov | ESTADO, data=int.beta.logit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ## Grafico na escala original p2 <- xyplot(logit(y) ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(preditos + lim.inf + lim.sup ~ cov | ESTADO, data=int.beta.logit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1,split=c(1,1,2,1),more=TRUE) plot(p2,split=c(2,1,2,1),more=FALSE) ########################################################################################## ## link Probit ########################################################################### ########################################################################################## ## Link probit inits <- coef(tra.probit) inits <- coef(pp.probit) beta.probit <- genericModel(formula = "y ~ cov + ESTADO", data = dados, inits = c(inits,1), family="beta", link = "probit") pp.probit <- profile(beta.probit,tol.newmin=1) plot(pp.probit) ## Plotando as bandas de confiança int.beta.probit <- cria.intervalo(beta.probit, formu= "~cov + ESTADO", data=newdados, nuisance=11) ## Plotando o grafico p1 <- xyplot(y ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(probit(preditos,inverse=TRUE) + probit(lim.inf,inverse=TRUE) + probit(lim.sup,inverse=TRUE) ~ cov | ESTADO, data=int.beta.probit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) p2 <- xyplot(probit(y) ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(preditos + lim.inf + lim.sup ~ cov | ESTADO, data=int.beta.probit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1,split=c(1,1,2,1),more=TRUE) plot(p2,split=c(2,1,2,1),more=FALSE) ########################################################################################## ## Link complemento complemento loglog ################################################### ########################################################################################## inits <- coef(tra.cloglog) inits <- coef(pp.cloglog) beta.cloglog <- genericModel(formula = "y ~ cov + ESTADO", data = dados, inits = c(inits,1), family="beta", link = "cloglog") pp.cloglog <- profile(beta.cloglog) plot(pp.cloglog) ## Plotando as bandas de confiança int.beta.cloglog <- cria.intervalo(beta.cloglog, formu= "~cov + ESTADO", data=newdados, nuisance=11) ## Plotando o grafico p1 <- xyplot(y ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(inv.cloglog(preditos) + inv.cloglog(lim.inf) + inv.cloglog(lim.sup) ~ cov | ESTADO, data=int.beta.cloglog, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) p2 <- xyplot(cloglog(y) ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(preditos + lim.inf + lim.sup ~ cov | ESTADO, data=int.beta.cloglog, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1,split=c(1,1,2,1),more=TRUE) plot(p2,split=c(2,1,2,1),more=FALSE) ########################################################################################## ## Complemento loglog #################################################################### ########################################################################################## inits <- coef(tra.loglog) inits <- coef(pp.loglog) beta.loglog <- genericModel(formula = "y ~ cov + ESTADO", data = dados, inits = c(inits,1), family="beta", link = "loglog") pp.loglog <- profile(beta.loglog) plot(pp.loglog) ## Plotando as bandas de confiança int.beta.loglog <- cria.intervalo(beta.loglog, formu= "~cov + ESTADO", data=newdados, nuisance=11) ## Plotando o grafico p1 <- xyplot(y ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(inv.loglog(preditos) + inv.loglog(lim.inf) + inv.loglog(lim.sup) ~ cov | ESTADO, data=int.beta.loglog, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) p2 <- xyplot(loglog(y) ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(preditos + lim.inf + lim.sup ~ cov | ESTADO, data=int.beta.loglog, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1,split=c(1,1,2,1),more=TRUE) plot(p2,split=c(2,1,2,1),more=FALSE) ########################################################################################## ## Link Cauchit ########################################################################## ########################################################################################## inits <- coef(tra.cauchit) inits <- coef(pp.cauchit) beta.cauchit <- genericModel(formula = "y ~ cov + ESTADO", data = dados, inits = c(inits,1), family="beta", link = "cauchit") pp.cauchit <- profile(beta.cauchit) plot(pp.cauchit) ## Plotando as bandas de confiança int.beta.cauchit <- cria.intervalo(beta.cauchit, formu= "~cov + ESTADO", data=newdados, nuisance=11) ## Plotando o grafico p1 <- xyplot(y ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(cauchit(preditos,inverse=TRUE) + cauchit(lim.inf,inverse=TRUE) + cauchit(lim.sup,inverse=TRUE) ~ cov | ESTADO, data=int.beta.cauchit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) p2 <- xyplot(cauchit(y) ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(preditos + lim.inf + lim.sup ~ cov | ESTADO, data=int.beta.cauchit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1,split=c(1,1,2,1),more=TRUE) plot(p2,split=c(2,1,2,1),more=FALSE) ## Link aranda inits <- coef(tra.logit) inits <- coef(pp.aranda) beta.aranda <- genericModel(formula = "y ~ cov + ESTADO", data = dados, inits = c(inits,1,1), family="beta", link = "aranda") pp.aranda <- profile(beta.aranda, tol.newmin=1) plot(pp.aranda) ## Plotando as bandas de confiança int.beta.aranda <- cria.intervalo(beta.aranda, formu= "~cov + ESTADO", data=newdados, nuisance=c(11,12)) ## Plotando o grafico p1 <- xyplot(y ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(inv.aranda(preditos,lambda=0.88) + inv.aranda(lim.inf,lambda=0.88) + inv.aranda(lim.sup,lambda=0.88) ~ cov | ESTADO, data=int.beta.cauchit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) p2 <- xyplot(aranda1(y,0.88) ~ cov | ESTADO,type="p", data=dados) + as.layer(xyplot(preditos + lim.inf + lim.sup ~ cov | ESTADO, data=int.beta.cauchit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1,split=c(1,1,2,1),more=TRUE) plot(p2,split=c(2,1,2,1),more=FALSE) ########################################################################################## ## Modelos Simplex ####################################################################### ########################################################################################## ## Link logit inits <- coef(tra.logit) inits <- coef(pp.simplex.logit) simplex.logit <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="simplex", link = "logit") pp.simplex.logit <- profile(simplex.logit) plot(pp.simplex.logit) ## Link probit inits <- coef(tra.probit) simplex.probit <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="simplex", link = "probit") pp.simplex.probit <- profile(simplex.probit) plot(pp.simplex.probit) ## Link complemento loglog inits <- coef(tra.cloglog) inits <- coef(pp.simplex.cloglog) simplex.cloglog <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="simplex", link = "cloglog") pp.simplex.cloglog <- profile(simplex.cloglog) plot(pp.simplex.cloglog) ## Link loglog inits <- coef(tra.loglog) simplex.loglog <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="simplex", link = "loglog") pp.simplex.loglog <- profile(simplex.loglog) plot(pp.simplex.loglog) ## Link Cauchit inits <- coef(tra.cauchit) simplex.cauchit <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="simplex", link = "cauchit") pp.simplex.cauchit <- profile(simplex.cauchit) plot(pp.simplex.cauchit) ## Link aranda inits <- coef(pp.simplex.aranda) simplex.aranda <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="simplex", link = "aranda") pp.simplex.aranda <- profile(simplex.aranda,tol.newmin=1) plot(pp.simplex.aranda) ########################################################################################## ## Modelos Kw ############################################################################ ########################################################################################## ## Link logit inits <- coef(tra.logit) inits <- coef(pp.kw.logit) kw.logit <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="Kw", link = "logit") pp.kw.logit <- profile(kw.logit) plot(pp.kw.logit) ## Link probit inits <- coef(tra.probit) inits <- coef(pp.kw.probit) kw.probit <- genericModel(formula = y ~ RENDA + PORTE + RAMO + ESTADO, data = dados, inits = inits, family="Kw", link = "probit") #VCOV <- vcov(kw.probit) pp.kw.probit <- profile(kw.probit,tol.newmin=1) plot(pp.kw.probit) ## Link complemento loglog inits <- coef(pp.kw.cloglog) kw.cloglog <- genericModel(formula = "y ~ RENDA + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="Kw", link = "cloglog") #VCOV <- vcov(kw.cloglog) pp.kw.cloglog <- profile(kw.cloglog) plot(pp.kw.cloglog) ## Link loglog inits <- coef(tra.loglog) inits <- coef(pp.kw.loglog) kw.loglog <- genericModel(formula = y ~ I(RENDA/100) + PORTE + RAMO + ESTADO, data = dados, inits = inits, family="Kw", link = "loglog") summary(kw.loglog) #VCOV <- vcov(kw.loglog) pp.kw.loglog <- profile(kw.loglog,tol.newmin=100,std.err = sqrt(diag(VCOV))) pp.kw.loglog <- profile(kw.loglog) plot(pp.kw.loglog) ## Link Cauchit inits <- coef(tra.cauchit) inits <- coef(pp.kw.cauchit) kw.cauchit <- genericModel(formula = "y ~ I(RENDA/100) + PORTE + RAMO + ESTADO", data = dados, inits = inits, family="Kw", link = "cauchit") summary(kw.cauchit) #VCOV <- vcov(kw.cauchit) pp.kw.cauchit <- profile(kw.cauchit,tol.newmin=3) plot(pp.kw.cauchit) ## Link aranda inits <- coef(tra.logit) inits <- coef(pp.kw.aranda) kw.aranda <- genericModel(formula = "y ~ I(RENDA/100) + PORTE + RAMO + ESTADO", data = dados, inits = c(inits,1,1), family="Kw", link = "aranda") pp.kw.aranda <- profile(kw.aranda) plot(pp.kw.aranda) ########################################################################################## ## Comparando os modelos vi log-verossimilhança ########################################## ########################################################################################## cbind(data.frame("LM" = logLik(mod.lm), "Lm.Logit" = ll.logit, "Lm.Probit" = ll.probit, "Lm.Cloglog" = ll.cloglog, "Lm.Loglog" = ll.loglog, "Lm.Cauchit" = ll.cauchit, "Lm.Aranda" = ll.aranda[1]), data.frame("B.Logit" = logLik(beta.logit), "B.Probit" = logLik(beta.probit), "B.Cloglog" = logLik(beta.cloglog), "B.Loglog" = logLik(beta.loglog),"B.Cauchit" = logLik(beta.cauchit), "B.Aranda" = logLik(beta.aranda)), data.frame("S.Logit" = logLik(simplex.logit), "S.Probit" = logLik(simplex.probit), "S.Cloglog" = logLik(simplex.cloglog), "S.Loglog" = logLik(simplex.loglog),"S.Cauchit" = logLik(simplex.cauchit), "S.Aranda" = logLik(simplex.aranda)), data.frame("Kw.Logit" = logLik(kw.logit), "Kw.Probit" = logLik(kw.probit), "Kw.Cloglog" = logLik(kw.cloglog), "Kw.Loglog" = logLik(kw.loglog),"Kw.Cauchit" = logLik(kw.cauchit), "Kw.Aranda" = logLik(kw.aranda)) ) ll.tra <- as.numeric(c(ll.logit, ll.probit, ll.cloglog, ll.loglog, ll.cauchit, ll.aranda[1])) ll.beta <- c(logLik(beta.logit),logLik(beta.probit),logLik(beta.cloglog),logLik(beta.loglog),logLik(beta.cauchit),logLik(beta.aranda)) ll.simplex <- c(logLik(simplex.logit),logLik(simplex.probit),logLik(simplex.cloglog),logLik(simplex.loglog),logLik(simplex.cauchit), logLik(simplex.aranda)) ll.kw <- c(logLik(kw.logit),logLik(kw.probit),logLik(kw.cloglog),logLik(kw.loglog),logLik(kw.cauchit),logLik(kw.aranda)) tab <- rbind(ll.tra,ll.beta,ll.simplex,ll.kw) tab colnames(tab) <- c("Logit", "Probit", "Cloglog", "Loglog","Cauchit", "Aranda") tab load("resultadosIQVT.RData") save.image("resultadosIQVT.RData")