# R code from vignette source 'cap06regnaolin.Rnw' # Encoding: UTF-8 #=========================================================================================== # code chunk number 1: cap06regnaolin.Rnw:3-4 #=========================================================================================== options(prompt=" ", continue=" ") #=========================================================================================== # code chunk number 2: cap06regnaolin.Rnw:11-71 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados de motivação # dia: dias após aplicação do nematicída # eclod: número de ovos eclodidos de nematóide lines <- " dia eclod 2 13.00 4 56.50 6 97.50 8 168.00 10 246.50 12 323.00 14 374.00 16 389.00 " da <- read.table(textConnection(lines), header=TRUE); closeAllConnections() str(da) plot(eclod~dia, da) #------------------------------------------------------------------------------------------ # ajuste de modelos lineares e não lineares new <- data.frame(dia=seq(0,30,l=100)) plot(eclod~dia, da, xlim=c(0,30), ylim=c(0,600)) #------------------------------------------------------------------------------------------ # modelo linear da reta m0 <- lm(eclod~dia, da) lines(predict(m0, newdata=new)~new$dia, col=1) #------------------------------------------------------------------------------------------ # modelo polinômio cúbico m1 <- lm(eclod~poly(dia, 3), da) lines(predict(m1, newdata=new)~new$dia, col=2) #------------------------------------------------------------------------------------------ # modelo não linear (logístico) m2 <- nls(eclod~SSlogis(dia, Asym, xmid, scal), data=da) lines(predict(m2, newdata=new)~new$dia, col=3) #------------------------------------------------------------------------------------------ # lineares: # * em termos de ajuste: é uma aproximação local e portanto a predição fora do intervalo # das covariáveis não deve ser feita # * em termos de interpretação: é um modelo empírico, relevância teórica, interpretação em # tem termos de taxa por unidade de incremento, nos polinômios interpretação é complicada # * em termos de inferência: estimadores de MQ0/verossimilhança possuem expressão analítica # inferências em pequenas amostras são exatas #------------------------------------------------------------------------------------------ # não lineares # * em termos de ajuste: é um modelo globa e portanto a predição fora do é justificada # * em termos de interpretação: é um modelo com obtido com suporte teórico, os parâmetros # possuem significado que agrega interpretação # * em termos de inferência: estimadores de MQ0/verossimilhança usam procedimentos # numéricos (gargalo do chute inicial), inferências em pequenas amostras são aproximadas e # dependem da parametrização do modelo e dos dados #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 3: cap06regnaolin.Rnw:76-94 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # as derivadas de y em relação a theta não são livres de theta # y = A*x/(B+x) : modelo michaelis mentem mm <- expression(A*x/(B+x)) # expressão do modelo D(mm, "A") # derivada simbólica de uma expressão D(mm, "B") #------------------------------------------------------------------------------------------ # y = A/(1+exp(B+C*x)) : modelo logístico logis <- expression(A/(1+exp(B+C*x))) D(logis, "A") D(logis, "B") D(logis, "C") mapply(D, name=c("A","B","C"), expr=logis) # calcula todas as derivadas de uma vez só #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 4: cap06regnaolin.Rnw:99-165 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # modelo michaelis mentem: f(x) = A*x/(B+x) # A: assintota (lim f(x) para x -> infinito) # B: tempo de meia vida (x para o qual f(x) = A/2) layout(1) A <- 10; B <- 3 curve(A*x/(B+x), 0, 50, ylim=c(0,10), col=2, lwd=3) abline(h=c(A, A/2), v=B, lty=3) #------------------------------------------------------------------------------------------ # se dividirmos numerador e denominador por B obtemos o ?monomolecular: f(x) = C*x/(1+x/B) # C: é a rezão entre assintota e tempo de meia vida pois C = A/B # B: é o tempo de meia vida C <- A/B curve(C*x/(1+x/B), add=TRUE, col=4, lwd=2, lty=2) # além de mudar intepretação (nenhuma intepretação cartesiana ou biológica para C), muda-se # as propriedades de estimação e inferência #------------------------------------------------------------------------------------------ # modelo logístico, f(x) = A/(1+exp((B-x)/C)) # A: assíntota (lim f(x) para x -> infinito) # B: tempo de meia vida (x para o qual f(x) = A/2) # C: ~escala, tamanho do intervalo de variação em x onde f(x) não é plana divido por 10 # o valor (0.09886*A/C) é a aproximadamente a taxa média nesse intervalo de x A <- 10; B <- 25; C <- 5 curve(A/(1+exp((B-x)/C)), 0, 50, col=2, lwd=3) abline(h=c(A, A/2), v=B, lty=3) rect(xleft=B-C*10/2, 0, B+C*10/2, A, density=5) abline(a=A*(-B+C*10/2)/(10*C), b=A/(10*C), col=3) xg <- seq(B-C*10/2, B+C*10/2,l=100) yg <- A/(1+exp((B-xg)/C)) mean(diff(yg)/diff(xg)) 0.09886*A/C #------------------------------------------------------------------------------------------ # modelo resposta platô: f(x) = A+B*x se x=x0 A <- 1; B <- 0.5; x0 <- 5 curve(A+B*x*(x=x0),0, 20, col=2, lwd=3) abline(h=c(A, A+B*x0), v=x0, lty=3) #------------------------------------------------------------------------------------------ # modelo de produção-competição (Bleasdale & Nelder, 1960): f(x) = x*(A+B*x)^(-1/C) # A: ? # B: ? # C: ? # o valor x = A/(B*(1/C-1)) é o ponto de máximo da função, para C<1 A <- 10; B <- 2; C <- 0.5 curve(x*(A+B*x)^(-1/C), 0, 50, col=2, lwd=3); abline(v=A/(B*(1/C-1))) C <- 1 curve(x*(A+B*x)^(-1/C), 0, 50, col=2, lwd=3) C <- 2 curve(x*(A+B*x)^(-1/C), 0, 50, col=2, lwd=3) #------------------------------------------------------------------------------------------ # modelo de curva de água (van Genuchten, 1980): f(x) = B+(A-B)/(1+(C*10^x)^n)^(1-1/n) A <- 0.7; B <- 0.3; C <- 1.3; D <- 1.6 curve(B+(A-B)/(1+(C*10^x)^D)^(1-1/D), -3, 4, col=2, lwd=3) abline(h=c(A,B), lty=3) par(new=TRUE) curve(eval(D(expression(B+(A-B)/(1+(C*10^x)^D)^(1-1/D)), "x")), -3, 4, yaxt="n", xaxt="n") axis(4) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 5: cap06regnaolin.Rnw:170-241 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # pacote que permite a construção de interfaces gráficas library(gWidgetsRGtk2) options("guiToolkit"="RGtk2") #------------------------------------------------------------------------------------------ # modelo michaelis mentem (reações químicas, liberação de nutrientes no solo) limits <- list(A=c(0,20), B=c(0,6)) plottest <- function(...){ curve(svalue(A)*x/(svalue(B)+x), 0, 15, ylim=c(0,10)) abline(v=svalue(B), h=svalue(A)/c(1,2), col=2) } #------------------------------------------------------------------------------------------ # função que constrói a janela gráfica com deslizadores #func <- ' w <- gwindow("Slider and spinbox example") tbl <- glayout(cont=w) for(i in 1:length(limits)){ tbl[i,1] <- paste("Slide to adjuste parameter", names(limits)[i]) tbl[i,2, expand=TRUE] <- (assign(names(limits)[i], gslider(from=limits[[i]][1], to=limits[[i]][2], by=diff(limits[[i]])/20, value=mean(limits[[i]]), container=tbl, handler=plottest))) } plottest() #'; writeLines(func, "func.R") #------------------------------------------------------------------------------------------ # modelo logístico (curva de crescimento) limits <- list(A=c(0,20), B=c(10,40), C=c(0.5,7)) plottest <- function(...){ AA <- svalue(A); BB <- svalue(B); CC <- svalue(C) curve(AA/(1+exp((BB-x)/CC)), 0, 50, ylim=c(0,15)) abline(h=c(AA, AA/2), v=BB, lty=3, col=2) rect(xleft=BB-CC*10/2, 0, BB+CC*10/2, AA, border="red") } source("func.R") #------------------------------------------------------------------------------------------ # modelo resposta platô (ensaios com fetilizante) limits <- list(A=c(0,2), B=c(0,2), x0=c(2,7)) plottest <- function(...){ AA <- svalue(A); BB <- svalue(B); x00 <- svalue(x0) curve(AA+BB*x*(x=x00), 0, 20, ylim=c(0,3)) abline(h=c(AA, AA+BB*x00), v=x00, col=2) } source("func.R") #------------------------------------------------------------------------------------------ # modelo de produção-competição (Bleasdale & Nelder, 1960) limits <- list(A=c(0,20), B=c(0,2), C=c(0,2)) plottest <- function(...){ AA <- svalue(A); BB <- svalue(B); CC <- svalue(C) curve(x*(AA+BB*x)^(-1/CC), 0, 50, ylim=c(0,0.7)) abline(v=AA/(BB*(1/CC-1)), col=2) } source("func.R") #------------------------------------------------------------------------------------------ # modelo de curva de água no solo (van Genuchten, 1980) limits <- list(A=c(0.5,1), B=c(0,0.5), C=c(0,5), D=c(1,6)) plottest <- function(...){ curve(svalue(B)+(svalue(A)-svalue(B))/(1+(svalue(C)*10^x)^svalue(D))^(1-1/svalue(D)), -3, 4, ylim=c(0,1)) } source("func.R") #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 6: cap06regnaolin.Rnw:246-314 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # como funciona o procedimento iterativo para estimar parâmetros? # exemplo com o modelo michaelis mentem e dados de mentirinha theta <- c(A=10, B=3) da <- data.frame(x=seq(1,20,2)) da$y <- theta["A"]*da$x/(theta["B"]+da$x)+rnorm(da$x,0,0.2) plot(y~x, da) #------------------------------------------------------------------------------------------ # sequência de estimativas até a convergência do procedimento de estimação # caminho pela superfície de mínimos quadrados sqe <- function(A, B, y, x){ hy <- (A*x)/(B+x); sum((y-hy)^2) } SQE <- Vectorize(sqe, c("A", "B")) A.grid <- seq(0,40,l=100) B.grid <- seq(0,20,l=100) sqe.surf <- outer(A.grid, B.grid, SQE, da$y, da$x) contour(A.grid, B.grid, sqe.surf, levels=(1:35)^2, xlab="A", ylab="B", col="gray70") #start.list <- locator(n=4); #start.list <- split(do.call(cbind, start.list), f=1:4) #start.list <- lapply(start.list, function(x){names(x) <- c("A","B"); x}) start.list <- list(s1=c(A=0.1,B=0.1), s2=c(A=40,B=20), s3=c(A=35,B=2.5), s4=c(A=18,B=18)) par(mfrow=c(2,2)) for(lis in 1:4){ contour(A.grid, B.grid, sqe.surf, levels=(seq(1,35,2))^2, xlab="A", ylab="B", col="gray70") sink("trace.txt") n0 <- nls(y~A*x/(B+x), data=da, start=start.list[[lis]], trace=TRUE) sink() trace <- read.table("trace.txt") for(i in seq(nrow(trace)-1)){ arrows(trace[i,"V3"], trace[i,"V4"], trace[i+1,"V3"], trace[i+1,"V4"], col=2, length=0.1) abline(v=trace[i+1,"V3"], h=trace[i+1,"V4"], col="orange", lty=3) Sys.sleep(1) print(c(i, trace[i+1,"V3"], trace[i+1,"V4"])) } } #------------------------------------------------------------------------------------------ # olhando a convergência pelo gráficos dos observados vs preditos for(lis in 1:4){ sink("trace.txt") n0 <- nls(y~A*x/(B+x), data=da, start=start.list[[lis]], trace=TRUE) sink() plot(y~x, da) trace <- read.table("trace.txt") for(i in seq(nrow(trace))){ curve(trace[i,"V3"]*x/(trace[i,"V4"]+x), add=TRUE, col=2) Sys.sleep(1) } } layout(1) #------------------------------------------------------------------------------------------ # curva ajustada plot(y~x, da) curve(coef(n0)["A"]*x/(coef(n0)["B"]+x), add=TRUE, col=2) #------------------------------------------------------------------------------------------ # estimativas dos parâmetros summary(n0) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 7: cap06regnaolin.Rnw:320-394 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # importando dados #dap <- read.table(file.choose(), header=TRUE) # no windows vai abrir janela de procura #dap <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/dap.txt", header=TRUE) dap <- read.table("../dados/dap.txt", header=TRUE) names(dap) <- c("d","h") str(dap) #------------------------------------------------------------------------------------------ # ordenando e tomando só os casos completos dap <- dap[order(dap$d),] dapcc <- dap[complete.cases(dap),] str(dapcc) #------------------------------------------------------------------------------------------ # análise gráfica exploratória dos dados plot(h~d, dapcc) #------------------------------------------------------------------------------------------ # análise gráfica do modelo candidato h = b0*(1-exp(b1*d))^b2 start <- list() # lista vazia irá armazenar os valores do último movimento de mouse limits <- list(b0=c(25,35), b1=c(0,0.5), b2=c(0.7, 1.3)) # amplitude de variação plottest <- function(...){ plot(h~d, dapcc) curve(svalue(b0)*(1-exp(-svalue(b1)*x))^svalue(b2), add=TRUE, col=2) start <<- list(b0=svalue(b0), b1=svalue(b1), b2=svalue(b2)) } source("func.R") start #------------------------------------------------------------------------------------------ # ajustar o modelo não linear (com os bons chutes do procedimento gráfico que só o R tem!) n0 <- nls(h~b0*(1-exp(-b1*d))^b2, data=dapcc, start=start, trace=TRUE) summary(n0) #------------------------------------------------------------------------------------------ # passando os chutes sem procedimento gráfico n0 <- nls(h~b0*(1-exp(-b1*d))^b2, data=list(A=35, B=0.1, C=1.3), start=start, trace=TRUE) summary(n0) #------------------------------------------------------------------------------------------ # ajustar o modelo não linear (com chutes sem noção, ver os tipos de mensagem de erro) n1 <- nls(h~b0*(1-exp(-b1*d))^b2, data=dapcc, start=list(b0=35, b1=0, b2=1.3), trace=TRUE) # b1=0: exp(0*x) = 1 constante n1 <- nls(h~b0*(1-exp(-b1*d))^b2, data=dapcc, start=list(b0=35, b1=-1, b2=1.3), trace=TRUE) # b1<0: problema na derivada numérica n1 <- nls(h~b0*(1-exp(-b1*d))^b2, data=dapcc, start=list(b0=35, b1=0.1, b2=-1), trace=TRUE) # vai para região singular #------------------------------------------------------------------------------------------ # verificação do ajuste plot(h~d, dapcc) lines(fitted(n0)~d, dapcc, col=2) #------------------------------------------------------------------------------------------ # não temos os gráficos de resíduos prontos para modelos não lineares, vamos contruí-los # extraindo valores r.cru <- residuals(n0) var(r.cru) r.pad <- residuals(n0, type="pearson") var(r.pad) fitd <- fitted(n0) #------------------------------------------------------------------------------------------ # fazendos os gráficos par(mfrow=c(1,3)) plot(r.cru~fitd) abline(h=0, lty=3) scatter.smooth(sqrt(abs(r.pad))~fitd) qqnorm(r.pad); qqline(r.pad, lty=2) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 8: cap06regnaolin.Rnw:399-455 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # intervalo de confiança para as estimativas confint.default(n0) # intervalo assintótico sempre é simétrico confint(n0) # observar os valores para entender o que o perfilhamento #------------------------------------------------------------------------------------------ # o intervalo de confiança perfilhado para um parâmetro prof <- profile(n0) prof$b0 layout(1) plot(prof$b0[,1]~prof$b0[,2][,1], type="l") abline(h=c(-1.96,0,1.96), v=c(29.84, 31.93, 36.96), col=2) #------------------------------------------------------------------------------------------ # o intervalo de confiânça de b2 contém o 1, será que preciso de b2? n1 <- nls(h~b0*(1-exp(-b1*d)), data=dapcc, start=list(b0=30, b1=0.1)) summary(n1) #------------------------------------------------------------------------------------------ # como ficou? layout(1) plot(h~d, dapcc) lines(fitted(n0)~d, dapcc, col=2) lines(fitted(n1)~d, dapcc, col=3) #------------------------------------------------------------------------------------------ # teste da razão de verossimilhança para H0: b2=1 anova(n1, n0) #------------------------------------------------------------------------------------------ # comparar o ajuste do modelo não linear com o linear escolhido na classe dos lineares -2*c(logLik(n1))+2*length(coef(n1)) #help(AIC, help_type="html") # AIC (k=2 parâmetros) # -2*c(logLik(m5))+2*3 # 966.5362 (k=3 parâmetros) -2*c(logLik(n1))+log(nrow(dapcc))*length(coef(n1)) # BIC (k=2 parâmetros) -2*c(logLik(m5))+length(coef(m5))*log(nrow(dapcc)) #------------------------------------------------------------------------------------------ # R2 em modelos não lineares (danger!, o modelo trivial deve ser um modelo aninhado) R2 <- function(nls.obj){ da <- eval(nls.obj$data) resp.name <- all.vars(summary(nls.obj)$formula)[1] names(da)[which(names(da)==resp.name)] <- "y" sqn <- deviance(nls.obj) sqe <- deviance(lm(y~1, da)) 1-(sqn/sqe) } R2(n0) R2(n1) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 9: cap06regnaolin.Rnw:460-480 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # extrai matrix de derivadas primeiras avaliada em theta, é como a matrix X dos lineares F <- attr(n1$m$fitted(), "gradient") F #------------------------------------------------------------------------------------------ # com essa matriz ajustamos um lm para ter gráficos de resíduos, !!remover intercepto!! m0 <- lm(h~-1+F, data=dapcc) # summary(), anova() não servem para nada #------------------------------------------------------------------------------------------ # gráfico de análise dos resíduos par(mfrow=c(2,2)) plot(m0) mtext("Análise de resíduos para modelo de regressão não linear", outer=TRUE, line=-2, cex=1.4) layout(1) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 10: cap06regnaolin.Rnw:485-547 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # curve de crescimento de aves para 2 sistemas de resfriamento de galpão frango <- expand.grid(dia=2:42, sistema=factor(c("A","B"))) frango$peso <- c( 80.18145, 89.98167, 132.14629, 192.04534, 167.68245, 191.45191, 220.74227, 212.98519, 230.82651, 346.32728, 391.14474, 407.79706, 441.54167, 499.63470, 575.36996, 603.35279, 678.09090, 763.96071, 787.66652, 921.68731, 959.13005, 1069.59008, 1150.70054, 1269.26359, 1313.35194, 1419.24574, 1532.63279, 1647.94630, 1722.91144, 1832.84384, 1921.09935, 1960.50372, 2062.17519, 2204.45014, 2258.73203, 2311.79432, 2466.26338, 2505.48039, 2521.81638, 2625.00725, 2728.60234, 201.41506, 240.71230, 289.29251, 215.56332, 294.79948, 297.17629, 346.07243, 358.03428, 393.36050, 388.47739, 477.51108, 420.89742, 490.44854, 605.53948, 629.18954, 659.28526, 713.87248, 773.69469, 887.45404, 943.04904, 970.29292, 980.20056, 1142.43274, 1197.28398, 1187.79456, 1243.54212, 1340.48431, 1453.78205, 1542.45519, 1596.08595, 1702.33500, 1801.46693, 1847.62131, 1860.69871, 2018.38835, 2046.97753, 2077.06034, 2236.60287, 2238.75234, 2302.30264, 2354.35641) #------------------------------------------------------------------------------------------ # análise gráfica exploratória require(lattice) xyplot(peso~dia, groups=sistema, data=frango, type=c("p","smooth")) xyplot(peso~dia|sistema, data=frango, type=c("p","smooth")) #------------------------------------------------------------------------------------------ # ajuste de curvas individuais com modelo logístico nA <- nls(peso~A/(1+exp(-(dia-d50)/S)), data=subset(frango, sistema=="A"), start=list(A=3000, d50=25, S=10)) summary(nA) #------------------------------------------------------------------------------------------ nB <- nls(peso~A/(1+exp(-(dia-d50)/S)), data=subset(frango, sistema=="B"), start=list(A=3000, d50=25, S=10)) summary(nB) #------------------------------------------------------------------------------------------ # fazer o ajuste das duas curvas num único nls(), estimativa do QMR é mais consistente nAB <- nls(peso~A[sistema]/(1+exp(-(dia-d50[sistema])/S[sistema])), data=frango, start=list( A=c(3200,3200), d50=c(28,30), S=c(8,10))) summary(nAB) #------------------------------------------------------------------------------------------ # fazer o gráfico dos valores ajustados/preditos new <- expand.grid(dia=0:70, sistema=factor(c("A","B"))) new$fit <- predict(nAB, newdata=new) #------------------------------------------------------------------------------------------ # gráfico layout(1) with(frango, plot(peso~dia, col=sistema, xlim=c(0,70), ylim=c(0,3200))) with(subset(new, sistema=="A"), lines(dia, fit)) with(subset(new, sistema=="B"), lines(dia, fit, col=2)) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 11: cap06regnaolin.Rnw:552-586 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # hipóteses: H0A: A1=A2, H0B: d501=d502, H0C: S1=S2 confint.default(nAB) # baseado em normalidade assintótica confint(nAB) # baseado em perfil de verossimilhança (modelo logístico bem "linear") #------------------------------------------------------------------------------------------ # testar H0A: ajustar um modelo em que A seja igual para os dois sistemas nAB2 <- nls(peso~A/(1+exp(-(dia-d50[sistema])/S[sistema])), data=frango, start=list( A=c(3200), d50=c(28,30), S=c(8,10))) summary(nAB2) #------------------------------------------------------------------------------------------ # empregar o teste da razão de verossimilhança para testar H0A: A1=A2 anova(nAB2, nAB) #------------------------------------------------------------------------------------------ # fazer o gráfico dos valores ajustados/preditos new <- expand.grid(dia=0:70, sistema=factor(c("A","B"))) new$fit <- predict(nAB2, newdata=new) #------------------------------------------------------------------------------------------ # gráfico layout(1) with(frango, plot(peso~dia, col=sistema, xlim=c(0,70), ylim=c(0,3200))) with(subset(new, sistema=="A"), lines(dia, fit)) with(subset(new, sistema=="B"), lines(dia, fit, col=2)) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 12: cap06regnaolin.Rnw:591-633 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # estimar a diferença de peso aos 42 entre os sistemas ("teste de médias") coef(nAB2) # estimativas vcov(nAB2) # matriz de covariância das estimativas #------------------------------------------------------------------------------------------ # "converter" o modelo para linear, permitido porque o inferência em não linear é linear F <- attr(nAB2$m$fitted(), "gradient") m0 <- lm(peso~-1+F, data=frango) coef(m0) # não faz muito sentido, mas é resultado de uma rotação round(vcov(nAB2),3) # modelo não linear original round(vcov(m0),3) # sua versão linear (mesma matriz) #------------------------------------------------------------------------------------------ # precisamos montar a matriz do contraste, ou seja, a F* para os valores x* model <- deriv3(~i*A/(1+exp((d501-dia)/S1))+(1-i)*A/(1+exp((d502-dia)/S2)), c("A","d501","d502","S1","S2","i"), function(A, d501, d502, S1, S2, i, dia){ NULL }) lA <- as.list(c(coef(nAB2), i=1, dia=42)) # sistema A aos 42 dias lB <- as.list(c(coef(nAB2), i=0, dia=42)) # sistema B aos 42 dias FA <- do.call(model, lA) # predito, F e H para sistema A aos 42 dias FB <- do.call(model, lB) # predito, F e H para sistema B aos 42 dias c(pesoA=FA[1], pesoB=FB[1]) # o peso dos animais aos 42 dias CC <- attr(FA, "gradient")[1:5]-attr(FB, "gradient")[1:5] # matriz do contraste #------------------------------------------------------------------------------------------ # estimativa do contraste CC%*%coef(m0) diff(predict(nAB2, newdata=data.frame(sistema=c("B","A"), dia=42))) #------------------------------------------------------------------------------------------ # usando a gmodels::estimable() require(gmodels) estimable(m0, rbind(H0D=CC), conf.int=0.95) #------------------------------------------------------------------------------------------ # e o R²? é para o modelo todo, e não um para cada curva R2(nAB2) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 13: cap06regnaolin.Rnw:638-694 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados de número de nematóides ecloditos em função dos dias e dose de nematícida nema <- expand.grid(dia=seq(2,16,2), dose=gl(4,1,labels=c(0,1,5,10))) nema$eclod <- c(13, 56.5, 97.5, 168, 246.5, 323, 374, 389, 7, 26, 64.5, 126, 207.5, 282, 334, 343, 5, 21.5, 45.5, 79, 118.5, 146, 167.5, 174.5, 3.25, 9.25, 12.5, 20.5, 32.25, 39.25, 40.25, 42.25) str(nema) xyplot(eclod~dia, groups=dose, data=nema, type="b", auto.key=TRUE) #------------------------------------------------------------------------------------------ # carrega o pacote nlme (do grupo dos recomendados) require(nlme) #------------------------------------------------------------------------------------------ # fazendo ajustes simultâneos, a SSlogis é uma selfStart nema <- groupedData(eclod~dia|dose, data=nema) # expressa a estrutura do dado nema$dose <- factor(nema$dose, ordered=FALSE) n0 <- nlsList(eclod~SSlogis(dia, Asym, xmid, scal), data=nema) summary(n0) plot(augPred(n0)) # faz o gráfico com preditos e observados coef(n0) #------------------------------------------------------------------------------------------ # na nlsList não tem como restrigir a estimação, então usamos a gnls # podemos passar o objeto nlsList ou montar sem atalhos gn0 <- gnls(eclod~SSlogis(dia, Asym, xmid, scal), data=nema, params=Asym+xmid+scal~dose, start=c(44,130,320,380, 7.3,1,1,1, 2.3,0,0,0)) anova(gn0, type="marginal") # é um teste de Wald, testa se o 3 desvios do intercepto=0 #------------------------------------------------------------------------------------------ # novos valores de dia para a predição de eclod new <- expand.grid(dia=seq(0,20,0.2), dose=levels(nema$dose)) new$eclod <- predict(gn0, newdata=new) xyplot(eclod~dia, groups=dose, data=new, type="l") #------------------------------------------------------------------------------------------ # incluir todos os resultados em um único gráfico tudo <- rbind(nema, new) tudo$tipo <- rep(c("obs","fit"), c(nrow(nema),nrow(new))) xyplot(eclod~dia|factor(dose), groups=tipo, data=tudo, distribute.type=TRUE, type=c("l","p","g"), main="O revisor que pede gráficos no Excel deve ir preso!", sub="Os gráficos do R são científicos!", xlab="Período após a aplicação dos nematicídas (dias)", ylab="Nematóides eclodidos", key=list(x=0.8, y=0.9, lines=list(lty=c(NULL,1), col=c("#0080ff","#ff00ff")), text=list(c("ajustado","observado"))), layout=c(4,1)#, scales=list(y="free") ) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 14: cap06regnaolin.Rnw:699-710 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # testar a hipótese H0: scal_i=scal para todo i, ou seja, scal comum gn1 <- gnls(eclod~SSlogis(dia, Asym, xmid, scal), data=nema, params=list(Asym+xmid~dose, scal~1), start=c(44,130,320,380, 7.3,1,1,1, 2.3)) anova(gn1, gn0) # testa a H0 por LRT, note p-valor diferente do anova(..., "marginal") #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 15: cap06regnaolin.Rnw:715-740 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados para a Curva de Retenção de Água, umidade (theta) e tensão (psi) cra <- data.frame(psi=c(0.01,1,2,4,6,8,10,33,60,100,500,1500,2787,4727,6840,7863, 9030,10000,10833,13070,17360,21960,26780,44860,69873,74623, 87287,104757,113817,147567,162723,245310,262217,298223), theta=c(0.779,0.554,0.468,0.406,0.373,0.36,0.344,0.309,0.298, 0.292,0.254,0.241,0.236,0.233,0.223,0.202,0.172,0.187,0.138, 0.098,0.07,0.058,0.052,0.036,0.029,0.0213,0.0178,0.0174, 0.0169,0.0137,0.0126,0.0109,0.0106,0.0053)) #------------------------------------------------------------------------------------------ # gráfico dos valores observados plot(theta~log10(psi), data=cra) #------------------------------------------------------------------------------------------ # função do modelo duplo van Genuchten, 7 parâmetros, cria a função do modelo f(x,theta) dvg <- function(x, ts, ti, tr, ae, at, ne, nt){ tr+(ti-tr)/((1+(at*10^x)^nt)^(1-1/nt))+(ts-ti)/((1+(ae*10^x)^ne)^(1-1/ne)) } dvg(log10(c(1,10,100,1000,10000)), # log 10 das tensões 0.7, 0.2, 0.05, 1.3, 0.0001, 1.5, 3.5) # valor dos parâmetros #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 16: cap06regnaolin.Rnw:745-814 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # obter os chutes dos 7 parâmetros via procedimento gráfico require(rpanel) # pacote com funções para manipulação gráfica start <- list() # cria uma lista vazia para receber os valores finais dvg.panel <- function(panel){ plot(panel$theta~log10(panel$psi)) curve(dvg(x, ts=panel$ts, ti=panel$ti, tr=panel$tr, ae=panel$ae, at=panel$at, ne=panel$ne, nt=panel$nt), add=TRUE) start <<- list(ts=panel$ts, ti=panel$ti, tr=panel$tr, ae=panel$ae, at=panel$at, ne=panel$ne, nt=panel$nt) panel } panel <- rp.control(theta=cra$theta, psi=cra$psi) rp.slider(panel, ts, 0.7, 0.9, initval=0.8, showvalue=TRUE, action=dvg.panel) rp.slider(panel, tr, 0, 0.1, initval=0.05, showvalue=TRUE, action=dvg.panel) rp.slider(panel, ti, 0.1, 0.5, initval=0.2, showvalue=TRUE, action=dvg.panel) rp.slider(panel, ae, 1, 3, initval=1.5, showvalue=TRUE, action=dvg.panel) rp.slider(panel, at, 0, 0.0001, initval=0.00005, showvalue=TRUE, action=dvg.panel) rp.slider(panel, ne, 1, 4, initval=1.5, showvalue=TRUE, action=dvg.panel) rp.slider(panel, nt, 2, 7, initval=4.5, showvalue=TRUE, action=dvg.panel) start #------------------------------------------------------------------------------------------ # usando a manipulate, que só está disponível para o RStudio require(manipulate) start <- list() # cria uma lista vazia para receber os valores finais manipulate({ plot(theta~log10(psi), data=cra) curve(dvg(x, ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt), add=TRUE) start <<- list(ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt) }, ts=slider(0.7, 0.9, initial=0.8), ti=slider(0.15, 0.25, initial=0.2), tr=slider(0, 0.10, initial=0.05), ae=slider(1.01, 3, initial=1.3), at=slider(0, 0.0001, initial=0.00005), ne=slider(1.01, 3, initial=1.65), nt=slider(1.8, 5, initial=4.3) ) start #------------------------------------------------------------------------------------------ # usando a gWidgetsRGtk2 library(gWidgetsRGtk2) options("guiToolkit"="RGtk2") limits <- list(ts=c(0.7,0.9), tr=c(0,0.1), ti=c(0.1,0.5), ae=c(1,3), at=c(0,0.0001), ne=c(1,4), nt=c(2,7)) plottest <- function(...){ plot(theta~log10(psi), data=cra) ts <- svalue(ts); tr <- svalue(tr); ti <- svalue(ti) ae <- svalue(ae); at <- svalue(at); ne <- svalue(ne); nt <- svalue(nt) curve(dvg(x, ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt), add=TRUE) start <<- list(ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt) } w <- gwindow("Slider and spinbox example") tbl <- glayout(cont=w) for(i in 1:length(limits)){ tbl[i,1] <- paste("Slide to adjuste parameter", names(limits)[i]) tbl[i,2, expand=TRUE] <- (assign(names(limits)[i], gslider(from=limits[[i]][1], to=limits[[i]][2], by=diff(limits[[i]])/20, value=mean(limits[[i]]), container=tbl, handler=plottest))) } plottest() start #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 17: cap06regnaolin.Rnw:819-853 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # start <- list(ts=0.772, ti=0.225, tr=0.011, ae=2.5861, at=0.0000788, ne=1.4637, nt=2.786) start # valores salvos do último movimento #------------------------------------------------------------------------------------------ # ajuste do modelo os dados usando os chutes do procedimento gráfico (muito fácil) n0 <- nls(theta~tr+(ti-tr)/((1+(at*psi)^nt)^(1-1/nt))+(ts-ti)/((1+(ae*psi)^ne)^(1-1/ne)), data=cra, start=start) summary(n0) # quadro de estimativas confint(n0) # intervalos de confiança perfilhados #------------------------------------------------------------------------------------------ # faz a diagnose dos resíduos qqnorm(residuals(n0)) # gráfico para normalidade plot(residuals(n0)~log10(cra$psi)) # gráfico para falta de ajuste plot(abs(residuals(n0))~fitted(n0)) # gráfico para homogeneidade de variância #------------------------------------------------------------------------------------------ # gráfico dos dados com a curva estimada lis <- c(list(x=NULL), as.list(coef(n0)), body(dvg)) plot(theta~log10(psi), data=cra, # faz o gráfico ylab=expression("Conteúdo de água no solo"~(theta~","~g~g^{-1})), # rótulo y xlab=expression("Tensão matricial"~(Psi~","~"kPa")), # rótulo x xaxt="n") tmp <- as.function(lis) curve(tmp, add=TRUE, col=2, lwd=1.5) # adiciona a curva axis(1, at=-2:5, label=as.character(10^(-2:5)), lwd.ticks=2) # escala log s <- log10(sort(sapply(1:9, function(x) x*10^(-3:6)))) axis(1, at=s, label=NA, lwd=0.5) # traços secundários abline(v=-2:6, h=seq(0,1,0.05), lty=3, col="gray50") # grade de referência #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 18: cap06regnaolin.Rnw:858-932 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados de postássio liberado em função do tempo klib <- data.frame(k=c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67, 91.28, 105.22, 72.74, 81.88, 97.62, 90.14, 89.88, 113.22, 90.91, 115.39, 112.63, 87.51, 104.69, 120.58, 114.32, 130.07, 117.65, 111.69, 128.54, 126.88, 127.00, 134.17, 149.66, 118.25, 132.67, 154.48, 129.11, 151.83, 147.66, 127.30), t=rep(c(15, 30, 45, 60, 75, 90, 120, 150, 180, 210, 240, 270), each=3)) #------------------------------------------------------------------------------------------ # criando função que representa o modelo exponencial reparametrizado, retorna f(x), F e H expo.der <- deriv3(~A*(1-exp(-log(2)*t/V))+D*t, c("A", "V", "D"), function(t, A, V, D) NULL) #------------------------------------------------------------------------------------------ # diagnose gráfica e primeiro chute plot(k~t, data=klib, xlab="Período de incubacão (dias)", ylab="Potássio liberado acumulado (mg/kg de solo)") A <- 90; V <- 20; D <- 0.2 curve(expo.der(x, A, V, D), add=TRUE, col=2) start <- list(A=A, V=V, D=D) #------------------------------------------------------------------------------------------ # ajustando o modelo aos dados a partir dos valores iniciais via gráfico n0 <- nls(k~expo.der(t, A, V, D), data=klib, start=start) summary(n0) confint(n0) #------------------------------------------------------------------------------------------ # valores preditos, gradiente e hessiano avaliado nos valores estimados str(n0) str(n0$m$fitted()) c(n0$m$fitted()) attr(n0$m$fitted(), "gradient") attr(n0$m$fitted(), "hessian") #------------------------------------------------------------------------------------------ # obtenção dos valores preditos pred <- data.frame(t=seq(0,300,l=100)) der <- do.call(expo.der, args=c(list(t=pred$t), as.list(coef(n0)))) F <- attr(der, "gradient") # gradiente avaliado no novo grid mais fino de t U <- chol(vcov(n0)) se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2))) # erro padrão ttabelado <- qt(c(.5, .025,.975), df=df.residual(n0)) # valor t para IC de 95% #------------------------------------------------------------------------------------------ # gráficos dos observados, preditos com IC, legenda e equações plot(k~t, data=klib, xlab="Período de incubacão (dias)", ylab="Potássio liberado acumulado (mg/kg de solo)", xlim=c(0,300), ylim=c(0,160)) matlines(pred$t, c(der)+outer(se, ttabelado), type="l", col=c(1,2,2), lty=c(1,2,2)) legend("bottomright", legend=c("valores observados", "valores preditos", "intervalo de confiança (95%)"), lty=c(NA,1,2), col=c(1,1,2), pch=c(1,NA,NA), bty="n") cf <- format(coef(n0), digits=3) text(par("usr")[1], par("usr")[4], adj=c(-0.05,1.5), label=substitute(hat(k)[total]==a%.%(1-e^{-ln(2)%.%t/v})+d%.%t, list(a=cf[1], v=cf[2], d=cf[3]))) abline(v=coef(n0)["V"], h=coef(n0)["A"], col="gray70") curve(expo.der(x, coef(n0)["A"], coef(n0)["V"], 0), add=TRUE, col=3) curve(expo.der(x, 0, 0, coef(n0)["D"]), add=TRUE, col=3) text(225, coef(n0)["A"], pos=3, label=substitute(hat(k)["fácil"]==a%.%(1-e^{-ln(2)%.%t/v}), list(a=cf[1], v=cf[2]))) text(173, 38, pos=3, srt=18, label=substitute(hat(k)["difícil"]==d%.%t, list(d=cf[3]))) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 19: cap06regnaolin.Rnw:937-1002 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # comparar 3 modelos: quociente, exponencial e exponencial reparametrizado (acima) # precisamos do F e H de todos eles expo.der <- deriv3(~A*(1-exp(-log(2)*t/V))+D*t, c("A", "V", "D"), function(t, A, V, D) NULL) exp.der <- deriv3(~A*(1-exp(-B*t))+D*t, c("A", "B", "D"), function(t, A, B, D) NULL) quo.der <- deriv3(~A*t/(V+t)+D*t, c("A", "V", "D"), function(t, A, V, D) NULL) #------------------------------------------------------------------------------------------ # ajustar os 3 modelos aos dados de liberação de potássio expo <- nls(k~expo.der(t, A, V, D), data=klib, start=list(A=90, V=20, D=0.2)) exp0 <- nls(k~exp.der(t, A, B, D), data=klib, start=list(A=A, B=0.05, D=D)) quo0 <- nls(k~quo.der(t, A, V, D), data=klib, start=list(A=90, V=20, D=0.2)) #------------------------------------------------------------------------------------------ # valores preditos por cada modelo pred <- data.frame(t=seq(0,300,l=100)) expop <- do.call(expo.der, args=c(list(t=pred$t), as.list(coef(expo)))) exp0p <- do.call(exp.der, args=c(list(t=pred$t), as.list(coef(exp0)))) quo0p <- do.call(quo.der, args=c(list(t=pred$t), as.list(coef(quo0)))) plot(k~t, data=klib, xlim=c(0,300), ylim=c(0,160)) matlines(pred$t, cbind(c(expop),c(exp0p),c(quo0p)), type="l", col=c(1,2,3), lty=c(1,2,3)) #------------------------------------------------------------------------------------------ # curvatura média de Bates & Watts, valor de corte é 0.3*sqrt(qf(19/20, p, n-p)) require(MASS) rms.curv(expo) rms.curv(exp0) # mesma curvatura intrínseca, maior curvatura devido à parametrização rms.curv(quo0) # menor curvatura intrínseca, intermediária devido à parametrização 0.3*sqrt(qf(19/20, length(coef(expo)), df.residual(expo))) #------------------------------------------------------------------------------------------ # vício de Box (função que eu desenvolvi na minha Dissertação de Mestrado) browseURL(URLencode("http://www.leg.ufpr.br/%7Ewalmes/docs/Walmes%20Marques%20Zeviani%20-%20Dissertacao.pdf")) biasbox <- function(nls.obj){ theta <- summary(nls.obj)$coef[,1] sd.theta <- summary(nls.obj)$coef[,2] F <- attr(nls.obj$m$fitted(), "gradient") H <- attr(nls.obj$m$fitted(), "hessian") sig <- summary(nls.obj)$sigma n <- dim(F)[1] FlFi <- t(F)%*%F d <- -(sig^2/2)*sapply(1:n, function(x){ sum(diag(solve(FlFi)%*%H[x, , ]))}) bias <- as.vector(solve(FlFi)%*%t(F)%*%d) names(bias) <- names(coef(nls.obj)) bias.sd <- 100*bias/sd.theta bias.th <- 100*bias/theta return(list("viés bruto"=bias, "%viés(theta)"=bias.th, "%viés(sd(theta))"=bias.sd)) } #------------------------------------------------------------------------------------------ # calcula o vício de box para as estimativas dos parâmetros biasbox(expo) biasbox(exp0) biasbox(quo0) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 20: cap06regnaolin.Rnw:1007-1040 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # é apenas uma outra forma de ajustar modelo, vantagem: obter contornos de verossimilhança fo <- function(th, x, y){ ## th: vetor numérico de parâmetros do modelo ## x: vetor da covariável ## y: vetor da resposta mu <- th[1]*(1-exp(-log(2)*x/th[2]))+th[3]*x sqr <- sum((y-mu)^2) s2 <- sqr/(length(y)) ll <- sum(dnorm(y, mu, sqrt(s2), log=TRUE)) ## retorna a log-verossimilhança return(ll) } #------------------------------------------------------------------------------------------ # avaliando a ll e a sqr com a função fobj() fo(coef(expo), x=klib$t, y=klib$k); logLik(expo) #------------------------------------------------------------------------------------------ # passando função para a optim() otimizar op <- optim(c(90, 20, 0.2), fo, x=klib$t, y=klib$k, method="BFGS", hessian=TRUE, control=list(fnscale=-1)) op$par # estimativas coef(expo) logLik(expo) # log-verossmilhança op$value solve(-op$hessian) # matriz de covariância das estimativas (denominador n) vcov(expo) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 21: cap06regnaolin.Rnw:1045-1077 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # o modelo A*(1-exp(-log(2)*x/V))+D*x é não linear em V apenas, se V for conhecido o # modelo assume forma linear A*z+D*x, em que z = (1-exp(-log(2)*x/V)) fo.conc <- function(th, x, y){ ## th: vetor numérico do parâmetro V ## x: vetor da covariável ## y: vetor da resposta y <- matrix(y) z <- (1-exp(-log(2)*x/th)) X <- model.matrix(~-1+z+x) ## calcula A e D via mínimos quadrados bt <- solve(crossprod(X), crossprod(X,y)) mu <- X%*%bt s2 <- crossprod(y-mu)/nrow(y) ll <- sum(dnorm(c(y), c(mu), sqrt(s2), log=TRUE)) ## retorna a verossimilhança attr(ll, "theta") <- c(A=c(bt[1,]), V=th, D=c(bt[2,])) return(ll) } fo.conc(coef(expo)[2], x=klib$t, y=klib$k) logLik(expo); coef(expo) #------------------------------------------------------------------------------------------ # estimação op <- optim(c(20), fo.conc, x=klib$t, y=klib$k, method="BFGS", hessian=TRUE, control=list(fnscale=-1)) op$par #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 22: cap06regnaolin.Rnw:1082-1115 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # fazendo o ajuste por busca em um universo discreto de valores (aproximar chute inicial) # com os procedimentos gráficos a nls2() perde importância em aproximar o chute inicial require(nls2) coef(expo) start.grid <- expand.grid(A=seq(40,120,length=10), V=seq(5,25,length=10), D=seq(0,0.5,length=10)) # cria um grid de 10*10*10=1000 valores formula(expo) n0.aprox <- nls2(formula(expo), data=klib, algorithm="grid-search", start=start.grid) # vai procurar o melhor dentro do grid summary(n0.aprox) coef(n0.aprox) # são valores iniciais que posso passar numa busca contínua com a nls() #------------------------------------------------------------------------------------------ # usando os valores ótimos retornados pelo "grid-search" na nls() n0 <- nls(formula(expo), data=klib, start=coef(n0.aprox)) summary(n0) # precisou de poucas interações cbind(coef(n0), coef(n0.aprox)) #------------------------------------------------------------------------------------------ # qual a real vantagem? passar as estimativas de vero concentrada para nls2() encontrar # os erros padrões, valor t, pvalor, IC... coef.op <- fo.conc(op$par, x=klib$t, y=klib$k) coef.op <- attr(coef.op, "theta") n0.op <- nls2(formula(expo), data=klib, start=coef.op, algorithm="grid-search") summary(n0.op) confint.default(n0.op) summary(n0) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 23: cap06regnaolin.Rnw:1120-1157 (eval = FALSE) #=========================================================================================== #------------------------------------------------------------------------------------------ # funções para análise de resíduos, gráficos para região de confiança e contornos de SQR require(nlstools) help(nlstools, help_type="html") #------------------------------------------------------------------------------------------ # modelo ajustado de classe nls class(expo) summary(expo) overview(expo) # summary com informação extra plotfit(expo) # observado vs ajustado points(klib$k, klib$k, pch=19, cex=0.3) points(klib$k, fitted(expo), pch=4, cex=0.3) #------------------------------------------------------------------------------------------ # gráficos resíduos~ajustados, autocorrelação, qqplot... par(mfrow=c(2,2)) plot(nlsResiduals(expo)) layout(1) #------------------------------------------------------------------------------------------ # os 6 gráficos par(mfrow=c(3,2)) for(i in 1:6) plot(nlsResiduals(expo), which=i) layout(1) #------------------------------------------------------------------------------------------ # contornos da soma de quadrado dos resíduos plot(nlsContourRSS(expo)) # demora um pouco, observar se os formatos são ellipticos plot(nlsContourRSS(expo), nlev=6, col=FALSE) #------------------------------------------------------------------------------------------ # IC baseado em proposição e aceitação dos valores internos à SQR (Beale1960) plot(nlsConfRegions(expo), bounds=TRUE) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 24: cap06regnaolin.Rnw:1162-1183 (eval = FALSE) #=========================================================================================== #------------------------------------------------------------------------------------------ # gráficos dos intervalos de confiânça baseados em perfil de verossimilhança require(MASS) par(mfrow=c(3,3)) plot(profile(expo)) plot(profile(exp0)) plot(profile(quo0)) layout(1) #------------------------------------------------------------------------------------------ # contornos de verossimilhança layout(1) plot(nlsContourRSS(expo, lseq=50), nlev=6, col=FALSE) x11() plot(nlsContourRSS(exp0, lseq=50), nlev=6, col=FALSE) x11() plot(nlsContourRSS(quo0, lseq=50), nlev=6, col=FALSE) # note: A e D são lineares -> contornos elípticos, os demais têm desvios de elipcidade #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 25: cap06regnaolin.Rnw:1188-1205 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # ajustando modelo parcialmente linear (procedimento iterativo apenas em V) n0 <- nls(k~cbind(A=1-exp(-log(2)*t/V), D=t), data=klib, algorithm="plinear", start=c(V=20)) summary(n0) #------------------------------------------------------------------------------------------ # resumo: # usando a nls() direto com chutes visuais # usando a optim() no modelo completo # usando a optim() no modelo concentrado/parcialmente linear # usando a nls() com chutes da nls2() # usando a nls2() com os "chutes" da optim() # usando a nls() com um modelo concentrado/parcialmente linear #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 26: cap06regnaolin.Rnw:1210-1279 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados de umidade em função do tempo (não é medida repetida) sec <- read.table("../dados/secagem.txt", header=TRUE, sep="\t") str(sec) #------------------------------------------------------------------------------------------ # exploração gráfica require(lattice) xyplot(umid~tempo|nome, data=sec, type=c("p","smooth")) #------------------------------------------------------------------------------------------ # os solos foram secos em estufa e adicionados 40% de umidade, será que o micro remove ela? xyplot(umrel~tempo|nome, data=sec, type=c("p","smooth","g")) #------------------------------------------------------------------------------------------ # fazer um intervalo de confiança para a umidade em 40 minutos, ajustar modelo logístico require(nlme) # permite um ajuste conjunto com mais facilidade sec <- groupedData(umrel~tempo|nome, data=sec) n0 <- nlsList(umrel~SSlogis(tempo, A, x0, S), data=sec) summary(n0) coef(n0) qqnorm(residuals(n0)); qqline(residuals(n0)) plot(residuals(n0)~fitted(n0)) plot(augPred(n0)) #------------------------------------------------------------------------------------------ # fazendo a predição, controlando intervalo, etc pred <- expand.grid(tempo=0:60, nome=levels(sec$nome)) pred$umrel <- predict(n0, newdata=pred) pred$tipo <- "predito" sec$tipo <- "observado" pred <- rbind(pred, sec[,c("nome","tempo","umrel","tipo")]) #------------------------------------------------------------------------------------------ # o gráfico xyplot(umrel~tempo|nome, groups=tipo, data=pred, distribute.type=TRUE, type=c("p","l")) xyplot(umrel~tempo, groups=nome, data=subset(pred, tipo=="predito"), type="l") #------------------------------------------------------------------------------------------ # H0: 40 minutos é suficiente para secar o solo lvadbw <- nls(umrel~SSlogis(tempo, A, x0, S), data=subset(sec, nome=="LVAd-Bw")) summary(lvadbw) #------------------------------------------------------------------------------------------ # obter a banda de confiança para a curva F <- attr(lvadbw$m$fitted(), "gradient") tc <- qt(0.975, df=df.residual(lvadbw)) vc <- vcov(lvadbw) se <- diag(sqrt(F%*%vc%*%t(F))) pred2 <- data.frame(fit=fitted(lvadbw)) pred2$lwr <- pred2$fit-tc*se pred2$upr <- pred2$fit+tc*se pred2$tempo <- sec$tempo[sec$nome=="LVAd-Bw"] pred2$umrel <- sec$umrel[sec$nome=="LVAd-Bw"] #------------------------------------------------------------------------------------------ # gráfico with(pred2, matplot(tempo, cbind(fit, lwr, upr), type="l", col=c(1,2,2), lty=c(1,2,2))) with(pred2, points(tempo, umrel)) abline(v=seq(0,45,2.5), h=seq(0,1.1,0.05), col="gray90", lty=2) abline(h=1, v=40) #------------------------------------------------------------------------------------------ # será que os alguns solos possuem o mesmo padrão de secamento? #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 27: cap06regnaolin.Rnw:1284-1347 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # a variância é função de x, f(x)? par(mfrow=c(1,2)) scatter.smooth(abs(residuals(n0))~fitted(n0), col=2) # decaimento scatter.smooth(abs(residuals(n0))~sec$tempo, col=2) # decaimento layout(1) #------------------------------------------------------------------------------------------ # ajustando modelo não linear modelando a variância help(varClasses, help_type="html") # modelos para a variância help(gnls, help_type="html") # função que possui argumento para ajuste da variância #------------------------------------------------------------------------------------------ # ajustando modelando a variância com estrutura varExp() coef(n0) sec$nome <- factor(sec$nome, ordered=FALSE) # por padrão é ordenado n1 <- gnls(umrel~SSlogis(tempo, A, x0, S), data=sec, params=A+x0+S~nome, start=list(1,0,0,0, 14,0,0,0, 7,0,0,0)) print(n1, corr=FALSE) #------------------------------------------------------------------------------------------ # ajustando varExp() n2 <- update(n1, weights=varExp(form=~tempo)) summary(n2) print(n2, corr=FALSE) anova(n2, n1) # testa a H0: expon=0 logLik(n2) #------------------------------------------------------------------------------------------ # ajustando varPower() n3 <- update(n1, weights=varPower(form=~tempo)) print(n3, corr=FALSE) anova(n3, n1) # testa a H0: power=0 logLik(n3) #------------------------------------------------------------------------------------------ # visualização gráfica das funções de variância aux <- tapply(residuals(n0), sec$tempo, var) plot(aux~unique(sec$tempo)) curve(0.1164^2*exp(2*(-0.03753)*x), add=TRUE, col=2) curve(0.2245^2*x^(-2*0.4998), add=TRUE, col=3) #------------------------------------------------------------------------------------------ # impacto nas estimativas e erros padrões dos parâmetros summary(n1)$tTable summary(n2)$tTable #------------------------------------------------------------------------------------------ # portanto tem algum impacto nos valores preditos e seus erros padrões pred <- expand.grid(tempo=0:60, nome=levels(sec$nome)) pred$umrel1 <- predict(n1, newdata=pred) pred$umrel2 <- predict(n2, newdata=pred) #------------------------------------------------------------------------------------------ # gráfico dos valores preditos require(latticeExtra) xyplot(umrel~tempo|nome, data=sec)+ as.layer(xyplot(umrel1~tempo|nome, data=pred, type="l", col=1))+ as.layer(xyplot(umrel2~tempo|nome, data=pred, type="l", col=2)) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 28: cap06regnaolin.Rnw:1352-1535 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados de crescimento de goiaba goi <- read.table("../dados/goiaba.txt", header=TRUE) str(goi) #------------------------------------------------------------------------------------------ # análise exploratória require(lattice) splom(goi[c("long","trans","peso","volume")], groups=goi$daa) # falsa noção de correlação xyplot(long~trans|daa, goi) # correlação xyplot(long~trans|daa, goi, scales="free", type=c("p","r")) #------------------------------------------------------------------------------------------ # mais análise exploratória require(reshape) aux <- melt(goi, id.vars="daa", measure.vars=c("long","trans","peso","volume"), variable="medidas") str(aux) xyplot(value~daa|medidas, data=aux, type=c("p","a"), scales="free") #------------------------------------------------------------------------------------------ # ajuste de um modelo para a variável peso em função de daa, sem modelar a variância # o modelo é um Gompertz reparametrizado para esse estudo n0 <- nls(peso~A-(A-D)*exp(-exp(C*(daa-B))), data=goi, start=c(A=200, C=0.09, B=105, D=7)) summary(n0) #------------------------------------------------------------------------------------------ # análises gráficas do ajuste require(nlstools) par(mfrow=c(2,2)) plot(profile(n0)) par(mfrow=c(3,2)) for(i in 1:6) plot(nlstools::nlsResiduals(n0), which=i) layout(1) #------------------------------------------------------------------------------------------ # modelo com inclusão da função de variâcia: var(e)=sigma^2*|ln{daa}|^{2*delta} require(nlme) n1 <- gnls(peso~A-(A-D)*exp(-exp(C*(daa-B))), data=goi, start=c(A=225, C=0.05, B=109, D=7), weights=varPower(0.1, form=~log(daa))) print(n1, corr=FALSE) #------------------------------------------------------------------------------------------ # teste para H0: power=0 anova(n1, n0) intervals(n1) # intervalos assintóticos para os parâmetros #------------------------------------------------------------------------------------------ # compara os valores das estimativas e erros padrões summary(n0)$coeff summary(n1)$tTable # erro padrão de A aumentou... #------------------------------------------------------------------------------------------ # qqplot e residuos vs daa qqmath(~residuals(n0, type="pearson")+residuals(n1, type="pearson")) xyplot(residuals(n0, type="pearson")+residuals(n1, type="pearson")~goi$daa, jitter.x=TRUE, pch=19) #------------------------------------------------------------------------------------------ # gráfico dos valores preditos require(latticeExtra) pred <- data.frame(daa=seq(10, 140, l=100)) pred$n0 <- predict(n0, newdata=pred) pred$n1 <- predict(n1, newdata=pred) xyplot(peso~daa, data=goi)+ as.layer(xyplot(n0~daa, data=pred, type="l", col=1))+ as.layer(xyplot(pred$n1~daa, data=pred, type="l", col=2)) #------------------------------------------------------------------------------------------ # bandas de confiança para a média invgomp <- deriv3(~A-(A-D)*exp(-exp(C*(daa-B))), c("A","C","B","D"), function(daa,A,B,C,D){NULL}) # bandas para o modelo homogêneo A0 <- do.call(invgomp, args=c(as.list(coef(n0)), list(daa=pred$daa))) F0 <- attr(A0, "gradient") t0 <- qt(c(f0=0.5, l0=0.025, u0=0.975), df.residual(n0)) U0 <- chol(vcov(n0)) S0 <- sqrt(apply(F0%*%t(U0), 1, function(x) sum(x^2))) Y0 <- apply(outer(S0, t0, "*"), 2, function(x) x+c(A0)) # bandas para o modelo heterogêneo A1 <- do.call(invgomp, args=c(as.list(coef(n1)), list(daa=pred$daa))) F1 <- attr(A1, "gradient") t1 <- qt(c(f1=0.5, l1=0.025, u1=0.975), df.residual(n0)) U1 <- chol(vcov(n1)) S1 <- sqrt(apply(F%*%t(U1), 1, function(x) sum(x^2))) Y1 <- apply(outer(S1, t1, "*"), 2, function(x) x+c(A1)) pred <- cbind(pred, Y0, Y1) str(pred) #------------------------------------------------------------------------------------------ # gráfico com as bandas de confiância para a média xyplot(peso~daa, data=goi, jitter.x=TRUE)+ as.layer(xyplot(f0+l0+u0~daa, data=pred, type="l", col=1, lwd=c(2,1,1)))+ as.layer(xyplot(f1+l1+u1~daa, data=pred, type="l", col=2, lwd=c(2,1,1))) #------------------------------------------------------------------------------------------ # funções gráficas para fazer região preenchida entre as bandas de confiança prepanel.biH <- function(x, y, ly, uy, subscripts, ...){ x <- as.numeric(x) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) list(ylim=range(uy, ly, finite=TRUE), xlim=range(x)) } panel.biH <- function(x, y, ly, uy, subscripts, colr, alpha=0.1, ...){ y <- as.numeric(y) x <- as.numeric(x) or <- order(x) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) panel.polygon(c(x[or], x[rev(or)]), c(ly[or], uy[rev(or)]), col=colr, alpha=alpha, border=NA) panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1) panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1) panel.xyplot(x, y, subscripts=subscripts, ...) } #------------------------------------------------------------------------------------------ # região entre as bandas de confiança é preenchida xyplot(peso~daa, data=goi, type="n")+ as.layer(with(pred, xyplot(f0~daa, type="l", col=4, lwd=2, ly=l0, uy=u0, colr=4, alpha=0.3, prepanel=prepanel.biH, panel=panel.biH)) )+ as.layer(with(pred, xyplot(f1~daa, type="l", col=3, lwd=2, ly=l1, uy=u1, colr=3, alpha=0.3, prepanel=prepanel.biH, panel=panel.biH)) )+ as.layer(xyplot(peso~daa, data=goi, pch=1, col=1, jitter.x=TRUE)) #------------------------------------------------------------------------------------------ # calculando os intervalos de confiança para a média amostral em cada daa mic <- aggregate(goi$peso, list(daa=goi$daa), function(x){ mean(x)+c(0,-1,1)*qt(0.975, df=length(x)-1)*sd(x)/sqrt(length(x)) }) #------------------------------------------------------------------------------------------ # funções gráficas para fazer intervalos (bigodes) de confiança prepanel.ciH <- function(x, y, ly, uy, subscripts=subscripts, ...){ x <- as.numeric(x) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) list(ylim = range(y, uy, ly, finite = TRUE)) } panel.ciH <- function(x, y, ly, uy, subscripts=subscripts, pch=pch, cola, ...){ panel.xyplot(x, y, subscripts=subscripts, pch=pch, ...) x <- as.numeric(x) y <- as.numeric(y) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) panel.arrows(x, ly, x, uy, col=cola, length=0.05, angle=90, code=3) } #------------------------------------------------------------------------------------------ # regiões de confiança e intervalos de confiança para a média amostral xyplot(peso~daa, data=goi, type="n")+ as.layer(with(pred, xyplot(f0~daa, type="l", col=4, lwd=2, ly=l0, uy=u0, colr=4, alpha=0.3, prepanel=prepanel.biH, panel=panel.biH)) )+ as.layer(with(pred, xyplot(f1~daa, type="l", col=3, lwd=2, ly=l1, uy=u1, colr=3, alpha=0.3, prepanel=prepanel.biH, panel=panel.biH)) )+ as.layer(with(mic, xyplot(x[,1]~daa, ly=x[,2], uy=x[,3], pch=5, col=1, cola=1, prepanel=prepanel.ciH, panel=panel.ciH)) ) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 29: cap06regnaolin.Rnw:1540-1650 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados do livro Regression Analysis and Its Applications - Bates & Watts # The Chloride data frame has 54 rows and 2 columns representing measurements of the chloride # ion concentration in blood cells suspended in a salt solution. # conc: a numeric vector giving the chloride ion concentration (%). # time: a numeric vector giving the time of the concentration measurement (min). conc <- c(17.3, 17.6, 17.9, 18.3, 18.5, 18.9, 19, 19.3, 19.8, 19.9, 20.2, 20.5, 20.6, 21.1, 21.5, 21.9, 22, 22.3, 22.6, 22.8, 23, 23.2, 23.4, 23.7, 24, 24.2, 24.5, 25, 25.4, 25.5, 25.9, 25.9, 26.3, 26.2, 26.5, 26.5, 26.6, 27, 27, 27, 27, 27.3, 27.8, 28.1, 28.1, 28.1, 28.4, 28.6, 29, 29.2, 29.3, 29.4, 29.4, 29.4) time <- c(2.45, 2.55, 2.65, 2.75, 2.85, 2.95, 3.05, 3.15, 3.25, 3.35, 3.45, 3.55, 3.65, 3.75, 3.85, 3.95, 4.05, 4.15, 4.25, 4.35, 4.45, 4.55, 4.65, 4.75, 4.85, 4.95, 5.05, 5.15, 5.25, 5.35, 5.45, 5.55, 5.65, 5.75, 5.85, 5.95, 6.05, 6.15, 6.25, 6.35, 6.45, 6.55, 6.65, 6.75, 6.85, 6.95, 7.05, 7.15, 7.25, 7.35, 7.45, 7.55, 7.65, 7.75) plot(conc~time) lines(spline(conc~time)) #------------------------------------------------------------------------------------------ # ajustando o modelo de regressão não linear fm1 <- nls(conc~Asym*(1-prop*exp(-exp(lrc)*time)), start=c(Asym=50, prop=0.6, lrc=log(0.25)), trace=TRUE) plot(resid(fm1)~fitted(fm1), type="b") # plot shows patterns in the residuals abline(h=0, lty=2, lwd=1) r <- residuals(fm1) plot(r[-1], r[-length(resid(fm1))]) # plot de e_{t} ~ e_{t-1} abline(h=0, lty=2, lwd=1) cor(r[-1], r[-length(resid(fm1))]) acf(r) pacf(r) #------------------------------------------------------------------------------------------ # veja os erros padrões summary(fm1)$coeff #------------------------------------------------------------------------------------------ # ajustando modelo declarando correlação require(nlme) fm2 <- gnls(conc~Asym*(1-prop*exp(-exp(lrc)*time)), start=c(Asym=50, prop=0.6, lrc=log(0.25)), correlation=corExp(0.9, form=~time)) fm2 summary(fm2) anova(fm2, fm1) #------------------------------------------------------------------------------------------ # praticamente o mesmo ajuste plot(conc~time) lines(fitted(fm1)~time) lines(fitted(fm2)~time) #------------------------------------------------------------------------------------------ # compare os erros padrões summary(fm1)$coeff # são menores porque correlação positiva induz menor variância summary(fm2)$tTable #------------------------------------------------------------------------------------------ # qual é a cara da matriz de correlação? tm <- seq(0,1,by=0.1) outer(tm, tm, function(x, y) round(exp(-abs(x-y)/0.2619), 2)) #------------------------------------------------------------------------------------------ # fazer o gráfico com IC time.g <- seq(min(time), max(time), l=200) mdl <- deriv3(~Asym*(1-prop*exp(-exp(lrc)*time)), c("Asym","prop","lrc"), function(time,Asym,prop,lrc){NULL}) A0 <- do.call(mdl, args=c(as.list(coef(fm1)), list(time=time.g))) F0 <- attr(A0, "gradient") t0 <- qt(c(f0=0.5, l0=0.025, u0=0.975), df=length(conc)-length(coef(fm1))) U0 <- chol(vcov(fm1)) S0 <- sqrt(apply(F0%*%t(U0), 1, function(x) sum(x^2))) Y0 <- apply(outer(S0, t0, "*"), 2, function(x) x+c(A0)) A1 <- do.call(mdl, args=c(as.list(coef(fm2)), list(time=time.g))) F1 <- attr(A1, "gradient") t1 <- qt(c(f1=0.5, l1=0.025, u1=0.975), df=length(conc)-length(coef(fm2))-1) U1 <- chol(vcov(fm2)) S1 <- sqrt(apply(F1%*%t(U1), 1, function(x) sum(x^2))) Y1 <- apply(outer(S1, t1, "*"), 2, function(x) x+c(A1)) #------------------------------------------------------------------------------------------ # gráfico com bandas spl <- spline(conc~time) #xyplot(conc~time)+ xyplot(conc~time, xlim=c(6,7), ylim=c(26.5,28.5))+ as.layer(xyplot(Y1[,1]~time.g, type="l", col=1, lwd=1, ly=Y1[,2], uy=Y1[,3], colr=2, alpha=0.1, prepanel=prepanel.biH, panel=panel.biH))+ as.layer(xyplot(Y0[,1]~time.g, type="l", col=4, lwd=1, ly=Y0[,2], uy=Y0[,3], colr=4, alpha=0.4, prepanel=prepanel.biH, panel=panel.biH))+ as.layer(xyplot(spl$y~spl$x, type="l", col=2, lwd=1)) #------------------------------------------------------------------------------------------ # predição com a geoR. (talvez excluir isso) require(geoR) da <- as.geodata(data.frame(x=time, y=0, r=r)) dag <- data.frame(x=time.g, y=0) pars <- c(c(fm2$sigma), exp(c(fm2$modelStruct$corStruct))) kg <- krige.conv(da, location=dag, krige=krige.control(cov.pars=pars)) xyplot(conc~time)+ as.layer(xyplot(Y1[,1]~time.g, type="l", col=1, lwd=1, ly=Y1[,2], uy=Y1[,3], colr=2, alpha=0.1, prepanel=prepanel.biH, panel=panel.biH))+ as.layer(xyplot(c(Y1[,1]+kg$predict)~time.g, type="l", col=2, lwd=1)) #------------------------------------------------------------------------------------------ #=========================================================================================== # code chunk number 30: cap06regnaolin.Rnw:1655-1736 (eval = FALSE) #=========================================================================================== #------------------------------------------------------------------------------------------ # dados de mineralização de nitrogênio # tempos de coleta, composto de lixo, esterco de galinha tp <- c(30, 45, 60, 75, 90, 105, 120, 150, 180, 210, 240, 270) cl <- c(24.11, 49, 71.01, 93.49, 103.97, 111.48, 118.72, 122.77, 132.42, 149.94, 157.4, 162.14) eg <- c(33.71, 63.91, 79.79, 104.1, 116.79, 131.17, 145.18, 154.46, 177.5, 185.14, 197.61, 198.95) #------------------------------------------------------------------------------------------ # exploração gráfica require(lattice) xyplot(cl+eg~tp, type="b") #------------------------------------------------------------------------------------------ # ajustando modelo para esterco de galinha - exponencial reparametrizado de 3 parâmetros require(nlme) n0 <- gnls(eg~A*(1-exp(-log(2)*(tp)/V))+D*tp, star=list(A=200, V=100, D=1)) summary(n0) # parâmetros não significativos!! plot(eg~tp); lines(fitted(n0)~tp) # estimação plausível pacf(residuals(n0)) # sem indícios de correlação acf(residuals(n0)) # sem indícios de correlação #------------------------------------------------------------------------------------------ # mas se fosse para declarar a correlação seria assim n1 <- update(n0, start=as.list(coef(n0)), correlation=corCAR1(0.5, form=~tp)) summary(n1) # ainda parâmetros não significativos!! plot(eg~tp); lines(fitted(n0)~tp) # estimação plausível lines(fitted(n1)~tp, col=2) # mesma predição anova(n0, n1) # correlação CAR1 não significativa #------------------------------------------------------------------------------------------ # os parâmetros são não significativos, mas o modelo é significativo comparado ao trivial anova(n0, lm(eg~1)) #------------------------------------------------------------------------------------------ # desconfiar de má especificação do modelo, deve faltar parâmetro tp.g <- seq(0, 300, by=2) plot(eg~tp, xlim=range(tp.g), ylim=c(0,max(eg))) lines(predict(n0, newdata=data.frame(tp=tp.g))~tp.g) # estimação plausível #------------------------------------------------------------------------------------------ # tentar outro modelo não linear, sem restrição y -> 0 quando x -> 0 n2 <- gnls(eg~A*(1-exp(-log(2)*(tp-t0)/V))+D*tp, star=list(A=200, V=100, D=1, t0=0)) summary(n2) # parâmetros significativos!! plot(eg~tp, xlim=range(tp.g), ylim=c(0,max(eg))) lines(predict(n0, newdata=data.frame(tp=tp.g))~tp.g) # estimação plausível lines(predict(n2, newdata=data.frame(tp=tp.g))~tp.g, col=2) # estimação BEM plausível #------------------------------------------------------------------------------------------ # veja o que a má especificação do modelo pode fazer! # a liberação de nitrogênio não começa no instante zero, demora um tempo para iniciar anova(n2, n0) #------------------------------------------------------------------------------------------ # plotar o gráfico com IC expt0 <- deriv3(~A*(1-exp(-log(2)*(tp-t0)/V))+D*tp, c("A","V","D","t0"), function(tp,A,V,D,t0){NULL}) A0 <- do.call(expt0, args=c(as.list(coef(n2)), list(tp=tp.g))) F0 <- attr(A0, "gradient") t0 <- qt(c(f0=0.5, l0=0.025, u0=0.975), df=length(eg)-length(coef(n2))) U0 <- chol(vcov(n2)) S0 <- sqrt(apply(F0%*%t(U0), 1, function(x) sum(x^2))) Y0 <- apply(outer(S0, t0, "*"), 2, function(x) x+c(A0)) #------------------------------------------------------------------------------------------ # gráfico com bandas xyplot(eg~tp, xlim=c(0,300), ylim=c(0,250))+ as.layer(xyplot(Y0[,1]~tp.g, type="l", col=2, lwd=2, ly=Y0[,2], uy=Y0[,3], colr=3, alpha=0.2, prepanel=prepanel.biH, panel=panel.biH)) #------------------------------------------------------------------------------------------ # qual a interpretação dos parâmetros summary(n2)$tTable trellis.focus("panel", 1, 1) panel.abline(v=coef(n2)["t0"], h=coef(n2)["A"], lty=2) panel.abline(a=0, b=coef(n2)["D"], lty=2) panel.curve(coef(n2)["A"]*(1-exp(-log(2)*(x-coef(n2)["t0"])/coef(n2)["V"])), lty=2) trellis.unfocus() #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 31: cap06regnaolin.Rnw:1741-1750 (eval = FALSE) #=========================================================================================== # http://www.ncfaculty.net/dogle/fishR/index.html # http://cran.r-project.org/web/packages/scapeMCMC/vignettes/dsc-vignette.pdf # http://cran.r-project.org/web/packages/gnm/index.html # http://cran.r-project.org/web/packages/NISTnls/index.html # http://cran.r-project.org/web/packages/nlreg/index.html # http://cran.r-project.org/web/packages/nlrwr/index.html # http://cran.r-project.org/web/packages/npde/index.html # http://cran.r-project.org/web/packages/NRAIA/index.html # http://cran.r-project.org/web/packages/HydroMe/index.html #=========================================================================================== # code chunk number 32: cap06regnaolin.Rnw:1755-1821 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados msa <- read.table("../dados/msacum.txt", header=TRUE) str(msa) names(msa) <- tolower(names(msa)) summary(msa) #------------------------------------------------------------------------------------------ # gráficos require(lattice) xyplot(mst~dae, groups=cult, data=msa, type="o") xyplot(mst~gde, groups=cult, data=msa, type="o") xyplot(mst~dsv, groups=cult, data=msa, type="o") #------------------------------------------------------------------------------------------ # ajustando curvas separadas, usando a nlsList e a selfStart SSlogis() require(nlme) msa <- groupedData(mst~dsv|cult, data=msa) n0 <- nlsList(mst~SSlogis(dsv,A,x0,S), data=msa) plot(residuals(n0)~fitted(n0)) qqnorm(residuals(n0)); qqline(residuals(n0)) summary(n0) plot(augPred(n0)) # gráfico com preditos e observados plot(intervals(n0)) # gráfico com intervalos assintóticos #------------------------------------------------------------------------------------------ # ajustar com a nls e obter o mesmo msa$cult <- factor(msa$cult, ordered=FALSE) n1 <- nls(mst~A[cult]/(1+exp(-(dsv-x0[cult])/S[cult])), data=msa, start=list(A=coef(n0)[,1],x0=coef(n0)[,2],S=coef(n0)[,3])) summary(n1) ci <- confint(n1) ci <- cbind(expand.grid(cult=levels(msa$cult), par=c("A","x0","S")), as.data.frame(ci)) names(ci)[3:4] <- c("L","U") ci #------------------------------------------------------------------------------------------ # gráfico com os IC perfilhados, são assimétricos require(latticeExtra) segplot(cult~L+U|par, scales=list(x="free"), data=ci, band=FALSE, center=coef(n1), layout=c(3,1)) #------------------------------------------------------------------------------------------ # ajustando com a retrição de que P e C tem o mesmo x0 msa$cultx0 <- factor(msa$cult) levels(msa$cultx0) <- c("D","CP","CP") n2 <- nls(mst~A[cult]/(1+exp(-(dsv-x0[cultx0])/S[cult])), data=msa, start=list(A=coef(n0)[,1],x0=c(0.87,0.97),S=coef(n0)[,3])) summary(n2) confint(n2) anova(n1, n2) #------------------------------------------------------------------------------------------ # fazer o gráfico com os preditos e observados range(msa$dsv) pred <- expand.grid(dsv=seq(0,2,l=30), cult=levels(msa$cult)) pred$cultx0 <- factor(pred$cult) levels(pred$cultx0) <- c("D","CP","CP") pred$y <- predict(n2, newdata=pred) xyplot(mst~dsv, groups=cult, data=msa)+ as.layer(xyplot(y~dsv, groups=cult, data=pred, type="l")) #------------------------------------------------------------------------------------------ .