[R-br] regressão polinomial: gráfico e valor máximo

Walmes Zeviani walmeszeviani em gmail.com
Terça Junho 21 12:42:13 BRT 2016


O ponto de máximo você aplica o cálculo. A expressão é simples.

Este conjunto de dados além de outros do Banzatto e Kronka e outras obras
(Pimental, Sonia Vieira, Zimmermann, mais de 15 obras e 350 datasets) então
sendo documentadas no pacote desenvolvimento pelo PET Estatística UFPR
(alerta de propaganda), o labestData (*lab*oratório de *est*atística):
https://github.com/pet-estatistica/labestData. Este dataset é chamado de
BanzattoQd7.2.1 no pacote (
https://github.com/pet-estatistica/labestData/blob/88f404857f39a1b047231d8a1aab8000fcb305b1/data-raw/BanzattoQd7.2.1.txt).
O meu CMR faz a leitura a partir do fonte dos dados no repositório. Se
instalar o pacote, terá à disposição datasets das obras nacionais de
estatística para usar em sala de aula e tutoriais.

rm(list = ls())

url <- paste0("https://raw.githubusercontent.com/pet-estatistica/",
              "labestData/88f404857f39a1b047231d8a1aab8000fcb305b1/",
              "data-raw/BanzattoQd7.2.1.txt")
BanzattoQd7.2.1 <- read.table(url, header = TRUE, sep = "\t")
str(BanzattoQd7.2.1)

# Polinômio ortogonal.
m0 <- lm(peso ~ poly(gesso, degree = 2), data = BanzattoQd7.2.1)

# Polinômio.
m0 <- lm(peso ~ poly(gesso, degree = 2, raw = TRUE), data = BanzattoQd7.2.1)
m0 <- lm(peso ~ gesso + I(gesso^2), data = BanzattoQd7.2.1)

coef(m0)

pred <- data.frame(gesso = seq(min(BanzattoQd7.2.1$gesso),
                               max(BanzattoQd7.2.1$gesso),
                               length.out = 30))
ci <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, ci)

ylim <- extendrange(c(BanzattoQd7.2.1$peso, ci))
# Do cálculo: f(x) = a * x^2 + b * x +c,
# então f'(x) = 0 é max/min = -b/(2 * a).
xmax <- -coef(m0)["gesso"]/(2 * coef(m0)["I(gesso^2)"])

plot(peso ~ gesso, data = BanzattoQd7.2.1)
with(pred, matlines(x = gesso, y = cbind(fit, lwr, upr),
                    col = c(1, 2, 2), lty = c(1, 2, 2)))
abline(v = xmax,
       h = predict(m0, newdata = list(gesso = xmax)),
       lty = 2, col = 3)


​À
​ disposição.
Walmes.​
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160621/8ccac149/attachment.html>


Mais detalhes sobre a lista de discussão R-br