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