#' ---
#' title: "LCE0602 - Estat?stica Experimental - Aula5"
#' author: ""
#' date: "09/2023"
#' ---
#' 
#' Ragazzi (1979) utilizou um experimento inteiramente 
#' casualizado com quatro repeti??es para estudar o efeito  
#' de 7 doses de gesso: 0,  50, 100, 150, 200, 250 e 300 kg/ha 
#' sobre diversas caracter?sticas do feijoeiro. 
#' Aqui foi avaliado  peso de 1000 sementes.



#' # Entrada dos dados para anÁlise

rm(list=ls(all=TRUE))
dose<- rep(seq(0,300,50),each=4)
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)
dados=data.frame(dose,peso)
head(dados)
(media=tapply(dados$peso,dados$dose,mean))


#' # Gr?fico de dispers?o
library(ggplot2)
ggplot(dados, aes(x = dose, y = peso)) +
  geom_point() +
  stat_summary(fun.y = mean, geom = "point", color = "red")

#' # Ajuste do modelo 
#' 
#' Assumindo-se que as pressuposi??es foram todas atendidas.
###############################
require(ExpDes.pt)
dic(dose,peso,quali=F)

#' # Gr?fico da curva ajustada
fun.2 <- function(x) {
  140.7839 + 0.2736*x -0.0008*x^2
}

ggplot(dados, aes(x = dose, y = peso)) +
  stat_summary(fun.y = mean, geom = "point", color = "red") +
  stat_function(fun = fun.2) +
  ylim(120,180)+
  xlab("dose de gesso (kg/ha)") +
  ylab("peso de 1000 sementes de feij?o (g)")




#' # Outra forma de ajustar o modelo

dados$Dose <- as.factor(dados$dose)

#' # ANOVA
mod.anova <- lm(peso ~ Dose, dados)
anova(mod.anova)

#' - Modelo de regress?o linear + Desvios de regress?o
mod.lin <- lm(peso ~ dose + Dose, dados)
anova(mod.lin)

#' - Modelo de regress?o quadr?tica + Desvios de regress?o
mod.quad <- lm(peso ~ dose + I(dose^2) + Dose, dados)
anova(mod.quad)

#' - Modelo selecionado
mod.aj <- lm(peso ~ dose + I(dose^2), dados)


#' ## Coeficientes estimados
coef(mod.aj)

#' ## Fun??o ajustada
fun.mod <- function(x){
  coef(mod.aj)[1] + coef(mod.aj)[2]*x + coef(mod.aj)[3]*x^2
}

ggplot(dados, aes(x = dose, y = peso)) +
  stat_summary(fun.y = mean, geom = "point", color = "red") +
  stat_function(fun = fun.mod) +
  ylim(120,180)+
  xlab("dose de gesso (kg/ha)") +
  ylab("peso de 1000 sementes de feij?o (g)")