Atividades e Materiais

Saudações. Está previsto para a disciplina de estudos dirigidos em estatística a abordagem do tópico Modelos Não Lineares. Para início de discussão os seguintes materiais devem ser lidos:

• Capítulo 8 Non-Linear and Smooth Regression do Modern Applied Statistics with S-plus (MASS);
• Apêndice do John Fox sobre regressão não linear com o R (Nonlinear Regression and Nonlinear Least Squares);
• Exemplo de ajuste de modelo não linear com o R da página do Prof Paulo Justiniano (Ajuste de modelos não lineares);
• Fazer o estudo dos dados abaixo. Propor um modelo não linear para os dados. Fazer o ajuste no R com as informações e rotinas vistas no materiais acima.

#-----------------------------------------------------------
# potássio liberado acumulado para o esterco de codorna
 
klib <- c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
          91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
          89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
          104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
          128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
          132.67, 154.48, 129.11, 151.83, 147.66, 127.30)
 
#-----------------------------------------------------------
# tempo em que foram feitas as coletas
 
tempo <- rep(c(15, 30, 45, 60, 75, 90,
               120, 150, 180, 210, 240, 270), each=3)
liber <- data.frame(tempo, klib)
 
require(lattice)
 
#-----------------------------------------------------------
# previa gráfica dos dados
 
xyplot(klib~tempo, data=liber,
       type=c("p", "smooth"), col=1,
       xlab="Período de incubação (dias)",
       ylab="Potássio liberado acumulado (mg/kg de solo)")
 
#-----------------------------------------------------------
# conjunto de dados com as médias das repetições e prévia
 
lmedio <- data.frame(tempo=unique(liber$tempo),
                     kmedio=tapply(liber$klib,
                       liber$tempo, mean))
 
xyplot(kmedio~tempo, data=lmedio,
       type=c("p", "smooth"), col=1,
       xlab="Período de incubação (dias)",
       ylab="Potássio liberado acumulado (mg/kg de solo)")
 
#-----------------------------------------------------------
# daqui em diante é com você...
 
#-----------------------------------------------------------
# exemplo de uso da optim()
 
x <- 1:9
A <- 5
B <- 1
y <- A*x/(B+x)+rnorm(x,0,0.1)
plot(y~x)
curve(A*x/(B+x), add=TRUE)
 
#-----------------------------------------------------------
# definição da função objetivo
 
fun.objetivo <- function(theta, y, x){
  sum((y-theta[1]*x/(theta[2]+x))^2)
}
 
#-----------------------------------------------------------
# escolha de valores iniciais
 
start <- c(3,0.5)
 
#-----------------------------------------------------------
# optimização da função objetivo
 
opt <- optim(start, fun.objetivo, y=y, x=x)
opt
 
curve(opt$par[1]*x/(opt$par[2]+x), add=TRUE, col=2)
 
#-----------------------------------------------------------
# usando outra função objetivo
 
fun.objetivo <- function(theta, y, x){
  n <- length(y)
  -(-n/2*log(2*pi)-n/2*log(theta[3])-
    sum((y-theta[1]*x/(theta[2]+x))^2/(2*theta[3])))
}
 
#-----------------------------------------------------------
# os chutes
 
start <- c(3,0.5,0.1)
 
#-----------------------------------------------------------
# optimização
 
opt <- optim(start, fun.objetivo, y=y, x=x)
opt
 
curve(opt$par[1]*x/(opt$par[2]+x), add=TRUE, col=3)
#-----------------------------------------------------------