########################################################################################## ## Modelo aditivo apenas por enquanto #################################################### ########################################################################################## ## Vamos usar os mesmos dados da OD OD <- read.table("OD.csv", header=TRUE, sep=",", dec=",") names(OD) <- c("Ensaio", "Rep", "Tipo", "Tempo", "OD") dados = subset(OD, Ensaio == c(1,2,4)) ## Lembrando dos dados require(lattice) require(latticeExtra) xyplot(OD ~ Tempo | factor(Ensaio), groups = Tipo, data=dados, type="p",layout=c(3,1)) ########################################################################################## ## Ajustando um modelo aditivo ########################################################### ########################################################################################## require(mgcv) ## Uma curva para todos os ensaios, repetições e tipo modelo <- gam(OD ~ s(Tempo), data=dados) plot(modelo) ## Uma curva para cada tipo modelo2 <- gam(OD ~ Tipo + s(Tempo, by=Tipo), data=dados) summary(modelo2) par(mfrow=c(2,1)) plot(modelo2) ## Uma curva para cada ensaio modelo3 <- gam(OD ~ factor(Ensaio) + s(Tempo, by=factor(Ensaio)), data=dados) summary(modelo3) par(mfrow=c(3,1)) plot(modelo3) ## Uma curva para cada ensaio e cada tipo dados$int <- interaction(dados$Ensaio, dados$Tipo) modelo4 <- gam(OD ~ factor(int) + s(Tempo, by=factor(int)), data=dados) summary(modelo4) par(mfrow=c(2,3)) plot(modelo4) ### Comparando os modelos para ver qual tem o melhor ajuste cbind(AIC(modelo),AIC(modelo2),AIC(modelo3),AIC(modelo4)) cbind(logLik(modelo),logLik(modelo2),logLik(modelo3),logLik(modelo4)) ## Predizendo com o modelo 4 Ensaio <- levels(factor(dados$Ensaio)) Tipo <- levels(factor(dados$Tipo)) Tempo <- seq(0,10, l=100) newdados <- expand.grid("Ensaio" = Ensaio, "Tipo" = Tipo, "Tempo" = Tempo) newdados$int <- interaction(newdados$Ensaio, newdados$Tipo) preditos <- predict(modelo4, newdados,type="response", se.fit=TRUE) preditos$sup <- preditos$fit + qnorm(0.975)*preditos$se.fit preditos$infe <- preditos$fit - qnorm(0.975)*preditos$se.fit xyplot(OD ~ Tempo | factor(Ensaio), groups = Tipo, data=dados, type="p",layout=c(3,1),auto.key=TRUE) + as.layer(xyplot(preditos$sup ~Tempo | factor(Ensaio), groups=Tipo, data=newdados, type="l")) + as.layer(xyplot(preditos$fit ~Tempo | factor(Ensaio), groups=Tipo, data=newdados, type="l")) + as.layer(xyplot(preditos$infe ~Tempo | factor(Ensaio), groups=Tipo, data=newdados, type="l")) xyplot(OD ~ Tempo | factor(Ensaio) + factor(Tipo), data=dados, type="p",layout=c(3,2),auto.key=TRUE) + as.layer(xyplot(preditos$sup ~Tempo | factor(Ensaio) + factor(Tipo), data=newdados, type="l")) + as.layer(xyplot(preditos$fit ~ Tempo | factor(Ensaio) + factor(Tipo), data=newdados, type="l")) + as.layer(xyplot(preditos$infe ~Tempo | factor(Ensaio) + factor(Tipo), data=newdados, type="l"))