########################################################################################## ## Comparando modelos de regressão para variáveis restrita ao intervalo (0,1) ############ ########################################################################################## ## Lendo a base de dados de severidade de doença dados <- read.table("dados.analise.txt", header=TRUE) ## Arrumando os ID dados$ID <- rep(1:6,times=length(unique(dados$dia))) ## Plotando um gráfico geral require(lattice) xyplot(inc ~ dia | cultivar, groups= ID,type="o", data=dados) ## Corrigindo observações na borda ## Corrigindo observações na borda 0 ou 1. dados$y <- (dados$inc/100*(length(dados$inc)-1) + 0.5)/length(dados$inc) ## Plotando novamente xyplot(y ~ dia | cultivar, groups= ID,type="o", data=dados) ## Pegando um subconjunto dadosU <- na.exclude(subset(dados, cultivar == "MARLI")) ## Plotando apenas esta cultivar xyplot(y ~ dia, groups=ID, type="o",data=dadosU) ########################################################################################## ## Ajustando alguns modelos ############################################################## ########################################################################################## require(bbmle) require(VGAM) require(latticeExtra) #require(rootSolve) ## Carregando as funções para os ajustes source("genericRegressao.R") source("functions.R") source("functionsKw.R") source("bandafunctions.R") ## family = beta newdados = data.frame(dia = seq(0,200,l=1000)) ## Link logit beta.logit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="beta", link = "logit") perf.logit <- profile(beta.logit) plot(perf.logit) ## Plotando as bandas de confiança int.beta.logit <- cria.intervalo(beta.logit, formu= "~ dia", data=newdados, nuisance=3, condicional=FALSE) ## Plotando o grafico p1 <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Logit", xlim=c(0,200)) + as.layer(xyplot(inv.logit(preditos) + inv.logit(lim.inf) + inv.logit(lim.sup) ~ dia, data=int.beta.logit, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Probit ########################################################################### ########################################################################################## beta.probit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="beta", link = "probit") perf.probit <- profile(beta.probit) plot(perf.probit) ## Plotando as bandas de confiança int.beta.probit <- cria.intervalo(beta.probit, formu= "~ dia", data=newdados, nuisance=3, condicional = FALSE) ## Plotando o grafico p2 <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Probit", xlim=c(0,200)) + as.layer(xyplot(probit(preditos,inverse=TRUE) + probit(lim.inf,inverse=TRUE) + probit(lim.sup,inverse=TRUE) ~ dia, data=int.beta.probit, type="l", col=c(2,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Complemento loglog #################################################################### ########################################################################################## ## link complemento log log beta.cloglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="beta", link = "cloglog") perf.cloglog <- profile(beta.cloglog) plot(perf.cloglog) ## Plotando as bandas de confiança int.beta.cloglog <- cria.intervalo(beta.cloglog, formu= "~ dia", data=newdados, nuisance=3, condicional=FALSE) ## Plotando o grafico p3 <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Comp. Log log",xlim=c(0,200)) + as.layer(xyplot(inv.cloglog(preditos) + inv.cloglog(lim.inf) + inv.cloglog(lim.sup) ~ dia, data=int.beta.cloglog, type="l", col=c(3,3,3), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## loglog ################################################################################ ########################################################################################## beta.loglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="beta", link = "loglog") perf.loglog <- profile(beta.loglog) plot(perf.loglog) ## Plotando as bandas de confiança int.beta.loglog <- cria.intervalo(beta.loglog, formu= "~ dia", data=newdados, nuisance=3, condicional=FALSE) ## Plotando o grafico p4 <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main= "Log log", xlim=c(0,200)) + as.layer(xyplot(inv.loglog(preditos) + inv.loglog(lim.inf) + inv.loglog(lim.sup) ~ dia, data=int.beta.loglog, type="l", col=c(4,4,4), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## link Cauchit ########################################################################## ########################################################################################## beta.cauchit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="beta", link = "cauchit") pp.cauchit <- profile(beta.cauchit) plot(pp.cauchit) ## Plotando as bandas de confiança int.beta.cauchy <- cria.intervalo(beta.cauchit, formu= "~ dia", data=newdados, nuisance=3,condicional=TRUE) ## Plotando o grafico p5 <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Cauchy", xlim=c(0,200)) + as.layer(xyplot(cauchit(preditos,inverse=TRUE) + cauchit(lim.inf,inverse=TRUE) + cauchit(lim.sup,inverse=TRUE) ~ dia, data=int.beta.cauchy, type="l", col=c(5,5,5), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Aranda Ordaz ##################################################################### ########################################################################################## ## Link Aranda beta.aranda <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1,1), family="beta", link = "aranda") pp.aranda <- profile(beta.aranda) plot(pp.aranda) ## Plotando as bandas de confiança int.beta.aranda <- cria.intervalo(beta.aranda, formu= "~ dia", data=newdados, nuisance=c(3,4),condicional=FALSE) ## Plotando o grafico p6 <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU,main="Aranda", xlim=c(0,200)) + as.layer(xyplot(inv.aranda(preditos,lambda=(coef(beta.aranda)[3])) + inv.aranda(lim.inf,lambda=(coef(beta.aranda)[3])) + inv.aranda(lim.sup,lambda=(coef(beta.aranda)[3])) ~ dia, data=int.beta.aranda, type="l", col=c(6,6,6), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Comparando todos os ajustes ########################################################### ########################################################################################## plot(p1,split=c(1,1,3,2),more=TRUE) plot(p2,split=c(2,1,3,2),more=TRUE) plot(p3,split=c(3,1,3,2),more=TRUE) plot(p4,split=c(1,2,3,2),more=TRUE) plot(p5,split=c(2,2,3,2),more=TRUE) plot(p6,split=c(3,2,3,2),more=FALSE) ########################################################################################## ## Comparando as verossimilhanca ######################################################### ########################################################################################## data.frame("Beta.Logit" = logLik(beta.logit), "Beta.Probit" = logLik(beta.probit), "cloglog" = logLik(beta.cloglog), "loglog" = logLik(beta.loglog),"cauchit" = logLik(beta.cauchit), "Aranda" = logLik(beta.aranda)) ########################################################################################## ## Ajustando com a simplex ############################################################### ########################################################################################## ########################################################################################## ## Simplex - Logit ####################################################################### ########################################################################################## inits <- c(0,0) inits <- coef(ss.logit) simplex.logit <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="simplex", link = "logit") ss.logit <- profile(simplex.logit) plot(ss.logit) ## Plotando as bandas de confiança int.simplex.logit <- cria.intervalo(simplex.logit, formu= "~ dia", data=newdados, nuisance=NULL,condicional=TRUE) ## Plotando o grafico p1.simplex <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Logit",xlim=c(0,200)) + as.layer(xyplot(inv.logit(preditos) + inv.logit(lim.inf) + inv.logit(lim.sup) ~ dia, data=int.simplex.logit, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Probit ########################################################################### ########################################################################################## inits <- c(0,0) inits <- coef(ss.probit) simplex.probit <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="simplex", link = "probit") ss.probit <- profile(simplex.probit) plot(ss.logit) ## Plotando as bandas de confiança int.simplex.probit <- cria.intervalo(simplex.probit, formu= "~ dia", data=newdados, nuisance=NULL, condicional=TRUE) ## Plotando o grafico p2.simplex <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Probit",xlim=c(0,200)) + as.layer(xyplot(probit(preditos,inverse=TRUE) + probit(lim.inf,inverse=TRUE) + probit(lim.sup,inverse=TRUE) ~ dia, data=int.simplex.probit, type="l", col=c(2,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## link Complemento LogLog ############################################################### ########################################################################################## inits <- c(0,0) inits <- coef(ss.cloglog) simplex.cloglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="simplex", link = "cloglog") ss.cloglog <- profile(simplex.cloglog) plot(ss.cloglog) ## Plotando as bandas de confiança int.simplex.cloglog <- cria.intervalo(simplex.cloglog, formu= "~ dia", data=newdados, nuisance=NULL, condicional=TRUE) ## Plotando o grafico p3.simplex <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Comp. Log log",xlim=c(0,200)) + as.layer(xyplot(inv.cloglog(preditos) + inv.cloglog(lim.inf) + inv.cloglog(lim.sup) ~ dia, data=int.simplex.cloglog, type="l", col=c(3,3,3), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Loglog ########################################################################### ########################################################################################## simplex.loglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0), family="simplex", link = "loglog") ss.loglog <- profile(simplex.loglog) plot(ss.loglog) ## Plotando as bandas de confiança int.simplex.loglog <- cria.intervalo(simplex.loglog, formu= "~ dia", data=newdados, nuisance=NULL,condicional=TRUE) ## Plotando o grafico p4.simplex <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main= "Log log",xlim=c(0,200)) + as.layer(xyplot(inv.loglog(preditos) + inv.loglog(lim.inf) + inv.loglog(lim.sup) ~ dia, data=int.simplex.loglog, type="l", col=c(4,4,4), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Cauchit ########################################################################## ########################################################################################## simplex.cauchit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0), family="simplex", link = "cauchit") ss.cauchit <- profile(simplex.cauchit) plot(ss.cauchit) ## Plotando as bandas de confiança int.simplex.cauchy <- cria.intervalo(simplex.cauchit, formu= "~ dia", data=newdados, nuisance=NULL,condicional=TRUE) ## Plotando o grafico p5.simplex <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Cauchy",xlim=c(0,200)) + as.layer(xyplot(cauchit(preditos,inverse=TRUE) + cauchit(lim.inf,inverse=TRUE) + cauchit(lim.sup,inverse=TRUE) ~ dia, data=int.simplex.cauchy, type="l", col=c(5,5,5), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## link Aranda ########################################################################### ########################################################################################## inits <- c(0,0,1) inits <- coef(ss.aranda) simplex.aranda <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="simplex", link = "aranda") ss.aranda <- profile(simplex.aranda) plot(ss.aranda) ## Plotando as bandas de confiança int.simplex.aranda <- cria.intervalo(simplex.aranda, formu= "~ dia", data=newdados, nuisance=3,condicional=FALSE) lambda = coef(simplex.aranda)[3] ## Plotando o grafico p6.simplex <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU,main="Aranda",xlim=c(0,200)) + as.layer(xyplot(inv.aranda(preditos,lambda=lambda) + inv.aranda(lim.inf,lambda=lambda) + inv.aranda(lim.sup,lambda=lambda) ~ dia, data=int.simplex.aranda, type="l", col=c(6,6,6), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Comparando os ajustes ################################################################# ########################################################################################## plot(p1.simplex,split=c(1,1,3,2),more=TRUE) plot(p2.simplex,split=c(2,1,3,2),more=TRUE) plot(p3.simplex,split=c(3,1,3,2),more=TRUE) plot(p4.simplex,split=c(1,2,3,2),more=TRUE) plot(p5.simplex,split=c(2,2,3,2),more=TRUE) plot(p6.simplex,split=c(3,2,3,2),more=FALSE) ########################################################################################## ## Comparando as verossimilhanca ######################################################### ########################################################################################## data.frame("Simplex.Logit" = logLik(simplex.logit), "simplex.Probit" = logLik(simplex.probit), "cloglog" = logLik(simplex.cloglog), "loglog" = logLik(simplex.loglog),"cauchit" = logLik(simplex.cauchit), "Aranda" = logLik(simplex.aranda)) ########################################################################################## ## family Kw ################################################################################ ################################################################################ ## Link logit kw.logit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="Kw", link = "logit") pkw.logit <- profile(kw.logit) plot(pkw.logit) ## Plotando as bandas de confiança int.kw.logit <- cria.intervalo(kw.logit, formu= "~ dia", data=newdados, nuisance=3,condicional=FALSE) ## Plotando o grafico p1.kw <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Logit",xlim=c(0,200)) + as.layer(xyplot(inv.logit(preditos) + inv.logit(lim.inf) + inv.logit(lim.sup) ~ dia, data=int.kw.logit, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link probit ########################################################################### ########################################################################################## kw.probit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="Kw", link = "probit") pkw.probit <- profile(kw.probit) plot(pkw.probit) ## Plotando as bandas de confiança int.kw.probit <- cria.intervalo(kw.probit, formu= "~ dia", data=newdados, nuisance=3, condicional=FALSE) ## Plotando o grafico p2.kw <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Probit",xlim=c(0,200)) + as.layer(xyplot(probit(preditos,inverse=TRUE) + probit(lim.inf,inverse=TRUE) + probit(lim.sup,inverse=TRUE) ~ dia, data=int.kw.probit, type="l", col=c(2,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link complemento log log ############################################################## ########################################################################################## inits <- coef(pkw.cloglog) inits <- c(-3,0.03,log(1)) kw.cloglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="Kw", link = "cloglog") pkw.cloglog <- profile(kw.cloglog) plot(pkw.cloglog) ## Plotando as bandas de confiança int.kw.cloglog <- cria.intervalo(kw.cloglog, formu= "~ dia", data=newdados, nuisance=3, condicional=FALSE) ## Plotando o grafico p3.kw <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Comp. Log log",xlim=c(0,200)) + as.layer(xyplot(inv.cloglog(preditos) + inv.cloglog(lim.inf) + inv.cloglog(lim.sup) ~ dia, data=int.kw.cloglog, type="l", col=c(3,3,3), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link loglog ########################################################################### ########################################################################################## kw.loglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="Kw", link = "loglog") pkw.loglog <- profile(kw.loglog) plot(pkw.loglog) ## Plotando as bandas de confiança int.kw.loglog <- cria.intervalo(kw.loglog, formu= "~ dia", data=newdados, nuisance=3,condicional=FALSE) ## Plotando o grafico p4.kw <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main= "Log log",xlim=c(0,200)) + as.layer(xyplot(inv.loglog(preditos) + inv.loglog(lim.inf) + inv.loglog(lim.sup) ~ dia, data=int.kw.loglog, type="l", col=c(4,4,4), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Cauchit ########################################################################## ########################################################################################## kw.cauchit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0,1), family="Kw", link = "cauchit") pkw.cauchit <- profile(kw.cauchit, tol.newmin=2) plot(pkw.cauchit) ## Plotando as bandas de confiança int.kw.cauchy <- cria.intervalo(kw.cauchit, formu= "~ dia", data=newdados, nuisance=3,condicional=FALSE) ## Plotando o grafico p5.kw <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Cauchy",xlim=c(0,200)) + as.layer(xyplot(cauchit(preditos,inverse=TRUE) + cauchit(lim.inf,inverse=TRUE) + cauchit(lim.sup,inverse=TRUE) ~ dia, data=int.kw.cauchy, type="l", col=c(5,5,5), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link aranda ########################################################################### ########################################################################################## inits <- c(0,0,1,1) inits <- coef(kw.aranda) kw.aranda <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="Kw", link = "aranda") pkw.aranda <- profile(kw.aranda) plot(pkw.aranda) ## Plotando as bandas de confiança int.kw.aranda <- cria.intervalo(kw.aranda, formu= "~ dia", data=newdados, nuisance=c(3,4),condicional=FALSE) ## Plotando o grafico lambda <- coef(kw.aranda)[3] p6.kw <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU,main="Aranda",xlim=c(0,200)) + as.layer(xyplot(inv.aranda(preditos,lambda=lambda) + inv.aranda(lim.inf,lambda=lambda) + inv.aranda(lim.sup,lambda=lambda) ~ dia, data=int.kw.aranda, type="l", col=c(6,6,6), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Comparando os ajustes ################################################################# ########################################################################################## plot(p1.kw,split=c(1,1,3,2),more=TRUE) plot(p2.kw,split=c(2,1,3,2),more=TRUE) plot(p3.kw,split=c(3,1,3,2),more=TRUE) plot(p4.kw,split=c(1,2,3,2),more=TRUE) plot(p5.kw,split=c(2,2,3,2),more=TRUE) plot(p6.kw,split=c(3,2,3,2),more=FALSE) ########################################################################################## ## Comparando as verossimilhanca ######################################################### ########################################################################################## data.frame("Kw.Logit" = logLik(kw.logit), "Kw.Probit" = logLik(kw.probit), "cloglog" = logLik(kw.cloglog), "loglog" = logLik(kw.loglog),"cauchit" = logLik(kw.cauchit), "Aranda" = logLik(kw.aranda)) ########################################################################################## ## family gaussian ####################################################################### ########################################################################################## ## link logit gs.logit <- genericModel(formula = "y ~ dia", data = dadosU, inits = c(0,0), family="gaussian", link = "logit") perf.gs.logit <- profile(gs.logit) plot(perf.gs.logit) ## Plotando as bandas de confiança int.gs.logit <- int.nao.linear(gs.logit, inv.link = inv.logitNL, formu = "~dia", newdados=newdados) ## Plotando o grafico p1.gs <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Logit") + as.layer(xyplot(preditos + lim.inf + lim.sup ~ dia, data=int.gs.logit, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Probit ################################################################################ ########################################################################################## inits <- coef(perf.gs.probit) inits <- c(0,0) gs.probit <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="gaussian", link = "probit") perf.gs.probit <- profile(gs.probit) plot(perf.gs.probit) ## Plotando as bandas de confiança int.gs.probit <- int.nao.linear(gs.probit, inv.link = inv.probitNL, formu = "~dia", newdados=newdados) ## Plotando o grafico p2.gs <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Probit") + as.layer(xyplot(preditos + lim.inf + lim.sup ~ dia, data=int.gs.probit, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Complemento loglog #################################################################### ########################################################################################## inits <- coef(perf.gs.cloglog) gs.cloglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="gaussian", link = "cloglog") perf.gs.cloglog <- profile(gs.cloglog) plot(perf.gs.cloglog) ## Plotando as bandas de confiança int.gs.cloglog <- int.nao.linear(gs.cloglog, inv.link = inv.cloglogNL, formu = "~dia", newdados=newdados) ## Plotando o grafico p3.gs <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "C loglog") + as.layer(xyplot(preditos + lim.inf + lim.sup ~ dia, data=int.gs.cloglog, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Loglog ################################################################################ ########################################################################################## inits <- coef(perf.gs.loglog) gs.loglog <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="gaussian", link = "loglog") perf.gs.loglog <- profile(gs.loglog) plot(perf.gs.loglog) ## Plotando as bandas de confiança int.gs.loglog <- int.nao.linear(gs.loglog, inv.link = inv.loglogNL, formu = "~dia", newdados=newdados) ## Plotando o grafico p4.gs <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Loglog") + as.layer(xyplot(preditos + lim.inf + lim.sup ~ dia, data=int.gs.loglog, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Cauchit ############################################################################### ########################################################################################## inits <- coef(perf.gs.cauchit) gs.cauchit <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="gaussian", link = "cauchit") perf.gs.cauchit <- profile(gs.cauchit) plot(perf.gs.cauchit) ## Plotando as bandas de confiança int.gs.cauchit <- int.nao.linear(gs.cauchit, inv.link = inv.cauchitNL, formu = "~dia", newdados=newdados) ## Plotando o grafico p5.gs <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Cauchit") + as.layer(xyplot(preditos + lim.inf + lim.sup ~ dia, data=int.gs.cauchit, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Aranda Ordaz ########################################################################## ########################################################################################## inits <- c(0,0,1) inits <- coef(perf.gs.aranda) gs.aranda <- genericModel(formula = "y ~ dia", data = dadosU, inits = inits, family="gaussian", link = "aranda") perf.gs.aranda <- profile(gs.aranda) plot(perf.gs.aranda) ## Plotando as bandas de confiança int.gs.aranda <- int.nao.linear(gs.aranda, inv.link = inv.arandaNL, formu = "~dia", newdados=newdados, lambda = coef(gs.aranda)[3]) ## Plotando o grafico p6.gs <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Aranda Ordaz") + as.layer(xyplot(preditos + lim.inf + lim.sup ~ dia, data=int.gs.aranda, type="l", col=c(1,1,1), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Comparando os ajustes ################################################################# ########################################################################################## plot(p1.gs,split=c(1,1,3,2),more=TRUE) plot(p2.gs,split=c(2,1,3,2),more=TRUE) plot(p3.gs,split=c(3,1,3,2),more=TRUE) plot(p4.gs,split=c(1,2,3,2),more=TRUE) plot(p5.gs,split=c(2,2,3,2),more=TRUE) plot(p6.gs,split=c(3,2,3,2),more=FALSE) ########################################################################################## ## Comparando as verossimilhanca ######################################################### ########################################################################################## data.frame("Gs.Logit" = logLik(gs.logit), "Gs.Probit" = logLik(gs.probit), "Gs.cloglog" = logLik(gs.cloglog), "Gs.loglog" = logLik(gs.loglog),"Gs.cauchit" = logLik(gs.cauchit), "Gs.Aranda" = logLik(gs.aranda)) ########################################################################################## ## modelos transformados ################################################################# ########################################################################################## tra.logit <- lm(logit(y) ~dia, data = dadosU) ll.logit <- sum(log(1/dadosU$y + 1/(1-dadosU$y))) + logLik(tra.logit) ## Bandas de confiança int.tra.logit <- int.tra(tra.logit, formu=" ~ dia", newdados=newdados, inv.link = inv.logit) ## Plotando o grafico p1.tra <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Logit") + as.layer(xyplot(predito + lim.inf + lim.sup ~ dia, data=int.tra.logit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) plot(p1.tra) ########################################################################################## ## Link Probit ########################################################################### ########################################################################################## tra.probit <- lm(probit(y) ~dia, data = dadosU) ll.probit <- sum(log(diag(rootSolve::gradient(qnorm, x = as.numeric(dadosU$y))))) + logLik(tra.probit) ## Bandas de confiança int.tra.probit <- int.tra(tra.probit, formu=" ~ dia", newdados=newdados, inv.link = inv.probit) ## Plotando o grafico p2.tra <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Probit") + as.layer(xyplot(predito + lim.inf + lim.sup ~ dia, data=int.tra.probit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## link Complemento Log Log ############################################################## ########################################################################################## tra.cloglog <- lm(cloglog(y) ~ dia, data = dadosU) ll.cloglog <- sum( log( - 1 / (log(1-dadosU$y)*(1-dadosU$y)))) + logLik(tra.cloglog) ## Bandas de confiança int.tra.cloglog <- int.tra(tra.cloglog, formu="~ dia", newdados=newdados, inv.link = inv.cloglog) ## Plotando o grafico p3.tra <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "C loglog") + as.layer(xyplot(predito + lim.inf + lim.sup ~ dia, data=int.tra.cloglog, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## link Log Log ########################################################################## ########################################################################################## tra.loglog <- lm(loglog(y) ~ dia, data = dadosU) ll.loglog <- sum(log( - 1/(dadosU$y*log(dadosU$y)))) + logLik(tra.loglog) ## Bandas de confiança int.tra.loglog <- int.tra(tra.loglog, formu=" ~ dia", newdados=newdados, inv.link = inv.loglog) ## Plotando o grafico p4.tra <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Loglog") + as.layer(xyplot(predito + lim.inf + lim.sup ~ dia, data=int.tra.loglog, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Cauchit ########################################################################## ########################################################################################## tra.cauchit <- lm(cauchit(y) ~ dia, data = dadosU) argum <- pi*(dadosU$y - 0.5) ll.cauchit <- sum(log(pi*(1/cos(argum))^2)) + logLik(tra.cauchit) ## Bandas de confiança int.tra.cauchit <- int.tra(tra.cauchit, formu=" ~ dia", newdados=newdados, inv.link = inv.cauchit) ## Plotando o grafico p5.tra <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Cauchit") + as.layer(xyplot(predito + lim.inf + lim.sup ~ dia, data=int.tra.cauchit, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Link Aranda Ordaz ##################################################################### ########################################################################################## mod.ordaz <- lm(y ~ dia, data=dadosU) ## Escolhendo o lambda ll.aranda <- aranda.ordaz(mod.ordaz, lambda = c(0.05,4)) tra.ordaz <- lm( aranda1(y, lambda = as.numeric(ll.aranda[2])) ~ dia, data=dadosU) ## Bandas de confiança int.tra.aranda <- int.tra(tra.ordaz, formu=" ~ dia", newdados=newdados, inv.link = inv.aranda, lambda = as.numeric(ll.aranda[2])) ## Plotando o grafico p6.tra <- xyplot(y ~ dia, groups=ID, type="o",data=dadosU, main = "Aranda Ordaz") + as.layer(xyplot(predito + lim.inf + lim.sup ~ dia, data=int.tra.aranda, type="l", col=c(1,2,2), lty = c(1,2,2), lwd = c(1,2,2))) ########################################################################################## ## Comparando os modelos transformados ################################################### ########################################################################################## plot(p1.tra,split=c(1,1,3,2),more=TRUE) plot(p2.tra,split=c(2,1,3,2),more=TRUE) plot(p3.tra,split=c(3,1,3,2),more=TRUE) plot(p4.tra,split=c(1,2,3,2),more=TRUE) plot(p5.tra,split=c(2,2,3,2),more=TRUE) plot(p6.tra,split=c(3,2,3,2),more=FALSE) ########################################################################################## ## Comparando as verossimilhanca ######################################################### ########################################################################################## data.frame("T.Logit" = ll.logit, "T.Probit" = ll.probit, "T.cloglog" = ll.cloglog, "T.loglog" = ll.loglog,"T.cauchit" = ll.cauchit, "T.Aranda" = ll.aranda[1]) ########################################################################################## ## Comparando os modelos ################################################################# ########################################################################################## ## Tabela de comparação ll.tra <- as.numeric(c(ll.logit, ll.probit, ll.cloglog, ll.loglog, ll.cauchit, ll.aranda[1])) ll.gs <- c(logLik(gs.logit), logLik(gs.probit), logLik(gs.cloglog), logLik(gs.loglog), logLik(gs.cauchit), logLik(gs.aranda)) 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.gs,ll.beta,ll.simplex,ll.kw) tab colnames(tab) <- c("Logit", "Probit", "Cloglog", "Loglog","Cauchit", "Aranda") tab ## 1 Cultivar save.image("modelosCORAL.RData") save(tab, file="CORAL.RData") ## 2 Cultivar save.image("modelosCHIRIPA.RData") save(tab, file="CHIRIPA.RData") ## 3 Cultivar save.image("modelosCHIMAR.RData") save(tab, file="CHIMAR.RData") ## 4 Cultivar save.image("modelosAURORA.RData") save(tab, file="AURORA.RData") ## 5 Cultivar save.image("modelosMACIEL.RData") save(tab, file="MACIEL.RData") ## 6 Cultivar save.image("modelosLEONENS.RData") save(tab, file="LEONENS.RData") ## 7 Cultivar save.image("modelosGRANADA.RData") save(tab, file="GRANADA.RData") ## 8 Cultivar save.image("modelosELDORA.RData") save(tab, file="ELDORA.RData") ## 9 Cultivar save.image("modelosVANGUA.RData") save(tab, file="VANGUA.RData") ## 10 Cultivar save.image("modelosPREMIER.RData") save(tab, file="PREMIER.RData") ## 11 Cultivar save.image("modelosMARLI.RData") save(tab, file="MARLI.RData") ########################################################################################## ## FIM DO CODIGO ######################################################################### ##########################################################################################