[R-br] regressão polinomial: gráfico e valor máximo (resolvido)
Maurício Lordêlo
mslordelo em gmail.com
Terça Junho 21 21:49:02 BRT 2016
Muito obrigado a todos.
Tudo esclarecido.
Maurício
Em 21 de junho de 2016 13:42, Éder Comunello <comunello.eder em gmail.com>
escreveu:
> Maurício e colegas, boa tarde!
>
> Mais uma opção de código para o problema.
> [image: Imagem inline 1]
>
> # <code r>
> # DIC 4 rep x 7 doses de gesso (trat): 0, 50, 100, 150, 200, 250, 300
> kg/ha
> # Peso de 1.000 sementes (peso, gramas)
> peso <- c(134.8, 139.7, 147.6, 132.3, 161.7, 157.7, 150.3, 144.7,
> 160.7, 172.7, 163.4, 161.3, 169.8, 168.2, 160.7, 161.0,
> 165.7, 160.0, 158.2, 151.0, 171.8, 157.3, 150.4, 160.4,
> 154.5, 160.4, 148.8, 154.0)
> trat <- rep(seq(0,300,50), each=4)
> dados <- data.frame(peso, trat=as.factor(trat))
>
> # Método dos Polinômios Ortogonais - Banzato & Kronka, p. 204
> contrasts(dados$trat) = contr.poly(levels(dados$trat));
> contrasts(dados$trat)
>
> fit1 = aov(peso ~ trat, dados)
> summary(fit1) # anova(fit1)
> summary(fit1, split = list(trat = list(LIN = 1, QUA = 2, CUB=3,
> DESVIOS=4:6)))
> # Df Sum Sq Mean Sq F value Pr(>F)
> # trat 6 1941.8 323.6 7.668 0.000188 ***
> # trat: LIN 1 423.2 423.2 10.026 0.004653 **
> # trat: QUA 1 1285.8 1285.8 30.465 1.78e-05 ***
> # trat: CUB 1 155.0 155.0 3.673 0.068997 .
> # trat: DESVIOS 3 77.8 25.9 0.614 0.613269
> # Residuals 21 886.3 42.2
>
> model.tables(fit1, "means")
> trat.m <- tapply(peso, trat, mean); trat.m
>
> # Uso do ajuste quadrático
> fit2 <- lm(peso ~ I(trat)+I(trat^2))
> summary(fit2)
> coef(fit2)
> eq <- paste0("Y = ", round(coef(fit2)[1], 4), " + ",
> round(coef(fit2)[2], 4), "X - ", abs(round(coef(fit2)[3],
> 6)), "X²"); eq
>
> # Coeficiente de determinação
> summary(fit2)$r.sq
> cor(fitted(fit2), peso)**2
> sumsq <- summary(fit1, split = list(trat = list(LIN = 1, QUA = 2, CUB=3,
> RES=4:6)))[[1]][2]
> rsq <- sum(sumsq[2:3,])/sumsq[1,]; rsq # Banzato & Kronka, p. 209
> eq.r <- paste0("R² = ", round(rsq, 4)); eq.r
>
> # Predições
> pred <- data.frame(trat = seq(0,300,1))
> pred$val <- predict(fit2, pred); head(pred)
>
> # Gráfico
> par(xpd=T)
> plot(c(0,300), c(135,170), type="n", xlab="Doses de gesso (kg/ha)",
> ylab="Peso de 1.000 sementes (g)", bty="n", xaxs="i", yaxs="i", las=1)
> points(as.numeric(names(trat.m)), trat.m, pch=18, cex=1)
> lines(pred$trat, pred$val, lwd=2, lty=1, col=1)
> text(300, 129, "X")
> text(-27, 170, "Y")
> text(150, 153, eq, adj=c(0.5,0.5))
> text(150, 151, eq.r, adj=c(0.5,0.5))
>
> # Max - Método 1
> ymax <- max(pred$val); ymax # 164.7042 gramas
> xmax <- pred[which.max(pred$val), "trat"]; xmax # com 175 Kg/ha de gesso
> segments(xmax, 135, xmax, ymax, lty=3, col=2)
> segments(0, ymax, xmax, ymax, lty=3, col=2)
> points(xmax, ymax, pch=20, col=2)
>
> # Max - Método 2
> # ax^2 + bx + c
> coefs <- data.frame(t(coef(fit2))); names(coefs) <- letters[3:1]; coefs
> with(coefs, -b/(2*a)) # 174.8403
> with(coefs, -(b^2 - 4*a*c)/(4*a)) # 164.7043
>
> # Max - Método 3
> fun <- function(x) 140.7839 + 0.2736*x - 0.000782*x**2
> fun(175)
> # optimize(fun, interval=c(0, 300), maximum=F)
> optimize(fun, interval=c(0, 300), maximum=T)
> # $maximum
> # [1] 174.9361
> #
> # $objective
> # [1] 164.7152
>
> # Max - Derivadas
>
> D1 <- D(expression(140.7839 + 0.2736*x - 0.000782*x**2), "x"); D1
> # 0.2736 - 0.000782 * (2 * x)
> # 0.2736 - 0.001564x
>
> D2 <- D(D1, "x"); D2
> # -(0.000782 * 2)
> # -0.001564
>
> # 0.2736 - 0.001564x = 0
> solve(0.001564, 0.2736) # Máximo
> # [1] 174.9361
> fun(174.9361) # 164.7152
>
> deriva <- deriv(~140.7839 + 0.2736*x - 0.000782*x**2,"x"); deriva
> x <- 175; eval(deriva)
> x <- 174.8403; eval(deriva)
> x <- 174.9361; eval(deriva)
> # </code>
>
>
> ================================================
> Éder Comunello
> Researcher at Brazilian Agricultural Research Corporation (Embrapa)
> DSc in Agricultural Systems Engineering (USP/Esalq)
> MSc in Environ. Sciences (UEM), Agronomist (UEM)
> ---
> Embrapa Agropecuária Oeste, Dourados, MS, Brazil |<O>|
> ================================================
> GEO, -22.2752, -54.8182, 408m
> UTC-04:00 / DST: UTC-03:00
>
>
>
>
> Em 20 de junho de 2016 19:34, Maurício Lordêlo <r-br em listas.c3sl.ufpr.br>
> escreveu:
>
>> Saudações a todos desta lista,
>> Estou tentando reproduzir um exemplo de regressão polinomial do livro do
>> Banzatto e Kronka.
>> Minhas dúvidas são:
>> 1. Como reproduzir exatamente o gráfico ilustrado no exemplo? O gráfico
>> que consegui fazer está parecido porém não está igual.
>> O gráfico apresentado no livro encontra-se neste link:
>> https://www.dropbox.com/s/qyv7ofwryegcu7m/figura.jpg?dl=0
>> 2. Como encontrar o valor de X que anula a derivada primeira? E depois
>> como encontrar o máximo da função?
>> Grato,
>> Maurício
>>
>> Segue o CMR:
>>
>> ################################################################################
>> #Experimento inteiramente casualizado com 4 repetições para estudar os
>> efeitos de
>> #7 doses de gesso (tratamentos):
>> #0, 50, 100, 150, 200, 250 e 300 kg/ha sobre diversas características do
>> feijoeiro
>> #Para a característica "peso de 1000 sementes", tem-se os resultados
>> obtidos em gramas
>> peso = c(134.8, 139.7, 147.6, 132.3,
>> 161.7, 157.7, 150.3, 144.7,
>> 160.7, 172.7, 163.4, 161.3,
>> 169.8, 168.2, 160.7, 161.0,
>> 165.7, 160.0, 158.2, 151.0,
>> 171.8, 157.3, 150.4, 160.4,
>> 154.5, 160.4, 148.8, 154.0)
>> tratamento = factor(rep(seq(0,300,50), each=4))
>> dados = data.frame(peso, tratamento)
>> attach(dados)
>>
>> niveis.tratamento = seq(0,300,50)
>> contrasts(dados$tratamento) = contr.poly(niveis.tratamento)
>> contrasts(dados$tratamento)
>> modelo.aov = aov(peso ~ tratamento, dados)
>> summary(modelo.aov)
>> summary(modelo.aov, split = list(tratamento = list(L = 1, Q = 2, C=3,
>> Dev=c(4,5,6))))
>>
>> plot(c(0,300), c(130,175), type="n", xlab="Doses de gesso (kg/ha)",
>> ylab="Peso (g) de 1000 sementes")
>> d = rep(seq(0,300,50), each=4)
>> med_t = rep(tapply(peso,tratamento,mean),each=4)
>> points(d,med_t)
>> x = seq(0,300,0.1)
>> lines(x, (predict(lm(peso ~ I(d)+I(d^2)+I(d^3)),data.frame(d=x))),lty=2,
>> col="blue")
>>
>>
>>
>> _______________________________________________
>> R-br mailing list
>> R-br em listas.c3sl.ufpr.br
>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>> código mínimo reproduzível.
>>
>
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160621/425d5509/attachment-0001.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: image.png
Tipo: image/png
Tamanho: 33792 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160621/425d5509/attachment-0001.png>
Mais detalhes sobre a lista de discussão R-br