#----------------------------------------------------------------------- # Ajuste de modelo não linear. library(lattice) library(latticeExtra) library(alr3) # Experimento de ganho de peso de perus com adição de metionina na # ração. plot(Gain ~ A, turk0) # Valores inicias para a curva. Foram determinados pela inspeção visual # do gráfico de dispersão. Int - intercepto, valor da resposta (eixo y) quando a # metionina é zero (eixo x), Amp - amplitude, tamanho do intervalo entro # o ponto no eixo y que é o intercepto e o ponto no eixo y que é # assíntota, Mei - "tempo de meia vida", corresponde ao valor no eixo x # onde a função vale Int + 0.5 * Amp. start <- list(Int = 625, Amp = 150, Mei = 0.08) plot(Gain ~ A, turk0) with(start, curve(Int + Amp * x/(Mei + x), add = TRUE, col = 2)) # Ajuste do modelo. n0 <- nls(Gain ~ Int + Amp * A/(Mei + A), data = turk0, start = start, trace = TRUE) # Gráfico da curva ajustada. plot(Gain ~ A, turk0) with(as.list(coef(n0)), curve(Int + Amp * x/(Mei + x), add = TRUE, col = 2)) # Ajuste de outro modelo. Trata-se de outra equação mas os parâmetros # tem o mesmo significado. n1 <- nls(Gain ~ Int + Amp * (1 - 2^(-A/Mei)), data = turk0, start = start, trace = TRUE) # Sobrepondo os dois ajustes. plot(Gain ~ A, turk0, xlim = c(0, 2), ylim = c(580, 850)) with(as.list(coef(n0)), curve(Int + Amp * x/(Mei + x), add = TRUE, col = 2)) with(as.list(coef(n1)), curve(Int + Amp * (1 - 2^(-x/Mei)), add = TRUE, col = 4)) legend("bottomright", legend = c("Michaelis-Menten Extendido", "Monomolecular reparametrizado"), lty = 1, col = c(2, 4), bty = "n") # Estimativas dos parâmentros. cbind(MME = coef(n0), MR = coef(n1)) # Resumo dos ajustes. Compare o valor do desvio-padrão residual. summary(n0) summary(n1) #----------------------------------------------------------------------- # Consumo de energia (KWH/dia) em função da temperatura (F). str(segreg) help(segreg, help_type = "html") # Passando de graus Fahrenheit para Celsius. segreg$cels <- (segreg$Temp * 5 - 160)/9 # Consumo médio mensal em função da temperatura média mensal. xyplot(C ~ cels, data = segreg) # Ajuste do modelo platô linear. # f(x) = Int + Taxa * (x - Queb) * (x >= Queb) # Int é o consumo de médio de energia antes do ponto limite ou ponto de # quebra Queb e Taxa é a taxa de aumento no consumo por unidade de # temperatura após o ponto de quebra. # Valores iniciais determinados visualmente. start <- list(Int = 75, Taxa = 1, Queb = 10) xyplot(C ~ cels, data = segreg) + layer(panel.curve(Int + Taxa * (x - Queb) * (x >= Queb)), data = start) # Ajuste do modelo. n0 <- nls(C ~ Int + Taxa * (cels - Queb) * (cels >= Queb), data = segreg, start = list(Int = 75, Taxa = 1, Queb = 10), trace = TRUE) # Modelo ajustado sobre os dados observados. xyplot(C ~ cels, data = segreg) + layer(panel.curve(Int + Taxa * (x - Queb) * (x >= Queb)), data = as.list(coef(n0))) # Estimativas dos parâmentros. summary(n0) # Intervalo de confiança (tipo Wald) para os parâmetros. confint.default(n0) # Modelo bi-linear onde as duas retas tem inclinação estimada. xyplot(C ~ cels, data = segreg) + layer(panel.curve((Int1 + Taxa1 * x) * (x <= Queb) + ((Int1 + Taxa1 * Queb) + Taxa2 * (x - Queb)) * (x > Queb)), data = list(Int1 = 75, Taxa1 = -0.05, Queb = 5, Taxa2 = 1)) # Ajuste do modelo bi-linear. n1 <- nls(C ~ (Int1 + Taxa1 * cels) * (cels <= Queb) + ((Int1 + Taxa1 * Queb) + Taxa2 * (cels - Queb)) * (cels > Queb), data = segreg, start = list(Int1 = 75, Taxa1 = -0.05, Queb = 5, Taxa2 = 1)) # Curva ajustada sobre os dados. xyplot(C ~ cels, data = segreg) + layer(panel.curve((Int1 + Taxa1 * x) * (x <= Queb) + ((Int1 + Taxa1 * Queb) + Taxa2 * (x - Queb)) * (x > Queb)), data = as.list(coef(n1))) # Estimativas dos parâmetros. summary(n1) # Teste os modelos encaixados: testa a hipótese de Taxa1 == 0. anova(n1, n0) #----------------------------------------------------------------------- # Tempos dos pitolos campeões (1º lugar) no círcuito de Monza. # Lendo tabela direto do site, requer o pacote XML. da <- XML::readHTMLTable("http://www.4mula1.ro/s/?cdx=9b0870e&pos=1&idt=44") da <- da[[1]] # Estrutura e aparência da tabela. str(da) head(da) # Anos das corridas. da$ano <- as.integer(sub("^(\\d{4})-.*$", "\\1", da$Date)) # Tempo do primeiro colocado na prova. library(lubridate) x <- hms(da$Time) da$tpro <- 60 * hour(x) + minute(x) + second(x)/60 # Cria anos a partir de 1950. da$ano <- da$ano - 1950 # Modelo bi-linear. xyplot(tpro ~ ano, data = da) + layer(panel.curve((Int1 + Taxa1 * x) * (x <= Queb) + ((Int1 + Taxa1 * Queb) + Taxa2 * (x - Queb)) * (x > Queb)), data = list(Int1 = 180, Taxa1 = -100/20, Queb = 20, Taxa2 = -5/30)) # Ajuste do modelo bi-linear. n1 <- nls(tpro ~ (Int1 + Taxa1 * ano) * (ano <= Queb) + ((Int1 + Taxa1 * Queb) + Taxa2 * (ano - Queb)) * (ano > Queb), data = da, start = list(Int1 = 180, Taxa1 = -100/20, Queb = 20, Taxa2 = -5/30)) # Curva ajustada. xyplot(tpro ~ ano, data = da) + layer(panel.curve((Int1 + Taxa1 * x) * (x <= Queb) + ((Int1 + Taxa1 * Queb) + Taxa2 * (x - Queb)) * (x > Queb)), data = as.list(coef(n1))) # Estimativas dos parâmetros. summary(n1) # Modelo logístico. start <- list(Ainf = 75, Asup = 180, x0 = 12, Sc = -8) xyplot(tpro ~ ano, data = da) + layer(panel.curve(Ainf + (Asup - Ainf)/(1 + exp(-(x - x0)/Sc))), data = start) # Ajuste do modelo. n2 <- nls(tpro ~ Ainf + (Asup - Ainf)/(1 + exp(-(ano - x0)/Sc)), data = da, start = start) summary(n2) # Modelo ajustado sobre os dados.. xyplot(tpro ~ ano, data = da) + layer(panel.curve(Ainf + (Asup - Ainf)/(1 + exp(-(x - x0)/Sc))), data = as.list(coef(n2))) # Desvio-padrão residual dos dois modelos. summary(n1)$sigma summary(n2)$sigma logLik(n1) logLik(n2) #-----------------------------------------------------------------------