Archive

Archive for the ‘regressão’ Category

Interface gráfica para ajuste de modelos de regressão não linear

Ajuste do modelo van Genuchten à curva de retenção de água do solo com interface gráfica desenvolvida com as ferramentas do pacote rpanel.

Modelos de regressão não linear são considerados quando alguma informação a priori existe sobre o fenômeno. Essa informação pode ser, por exemplo, que a curva seja sempre crescente, típico para curvas de crescimento/acúmulo. Em geral, esses modelos têm parâmetros com interpretação física/química/biológica e alguns parâmetros não têm, mas estão presentes para conferir flexibilidade.

Talvez um dos grandes gargalos no ajuste de modelos não lineares seja atribuição de valores iniciais para o método numérico de estimação. Para parâmetros com interpretação é possível obter valores ao ver gráficos de dispersão das variáveis ou pelo menos se tem idéia do intervalo em que o parâmetro se encontra. Para parâmetros sem interpretação a obtenção de valores iniciais é realmente trabalhosa.

Um procedimento muito útil para obter valores iniciais é fazer um grid de valores iniciais e plotar a curva corresponde para cada conjunto sobre o diagrama de dispersão dos dados. Você usa como valores iniciais os correspondentes à curva que melhor se aproxima do comportamento dos dados. Por outro lado, criar um grid de valores, sobrepor cada curva ao gráfico, escolher à curva que melhor se aproxima aos dados, fornecer os valores iniciais e estimar os parâmetros é algo de demanda tempo se o número de curvas à ajustar for grande. Seria de grande ajuda se esse processo pudesse ser automatizado ou facilitado, por exemplo, por meio de uma interface gráfica com botões de ação e controle de parâmetros por meio do mouse.

O pacote rpanel possui um conjunto de funções básicas que permitem construir interfaces gráficas para o usuário interagir com o R. O conjunto de funções compreende desde caixas de texto, caixas de seleção, deslizadores, até botões para disparo de alguma função. Usando essas ferramentas foi possível construir uma interface que permite selecionar o conjunto de dados, pré-ajustar a curva aos dados por meio de deslizadores e ajustar o modelo de regressão não linear. O código abaixo fornece as ferramentas para obter o que se vê nos gifs animados dessa matéria. Modificações devem ser feitas para outros dados/modelos. Minha intenção no futuro é aprimorar essas funções de forma que o usuário forneça o modelo e os limites dos intervalos para os parâmetros diretamente pela interface. Em caso de sucesso, isso pode virar um pacote ou complementar algum dos pacotes existentes dedicados à regressão não linear.

O primeiro exemplo usa gráficos do pacote graphics, básico do R. Os dados são de 10 curvas de retenção de água no solo onde ajustou-se o modelo van Genuchten (1980). O segundo conjunto de dados é de um experimento para avaliar o impacto da desfolha na produtividade do algodão. Para esse caso adotou-se o modelo potência e gráficos do pacote lattice. Clique na caixa abaixo para ver o código. Até à próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(rpanel)
require(lattice)
require(latticeExtra)
require(grid)

#-----------------------------------------------------------------------------
# dados de curva de retenção de água no solo

cra <- read.table("http://www.leg.ufpr.br/~walmes/data/cra_manejo.txt",
                  header=TRUE, sep="\t")
str(cra)
cra$tens[cra$tens==0] <- 0.1

#-----------------------------------------------------------------------------
# ver

cras <- subset(cra, condi=="LVA3,5")
cras <- with(cras, aggregate(cbind(umid),
                             list(posi=posi, tens=tens, prof=prof), mean))
cras$caso <- with(cras, interaction(posi, prof))

xyplot(umid~log10(tens)|posi, groups=prof, data=cras, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste do modelo van Genuchten de modo iterativo com rpanel

# nlsajust é excutada quando clicar no botão "Ajustar"
nlsajust <- function(panel){
  ## ajuste do modelo não linear
  n0 <- try(nls(umid~tr+(ts-tr)/(1+(a*tens)^n)^(1-1/n), # modelo
                data=da,                                # dados
                start=start))                           # valores iniciais
  ## em caso de não convergência imprime mensagem de erro
  if(class(n0)=="try-error"){
    par(usr=c(0, 1, 0, 1))
    text(0.5, 0.5, "Não convergiu!\nAproxime mais.", col="red", cex=2)
  } else {
    ## coloca à curva ajusta sobre os pontos
    with(as.list(coef(n0)), curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n),
                                  add=TRUE, col=2))
    ## salva o ajuste numa lista
    aju[[i]] <<- n0
  }
  panel
}

vg <- function(panel){
  ## seleciona o subconjunto dos dados
  i <<- panel$caso
  da <<- subset(cras, caso==i)
  ## lista com valores iniciais vindos dos deslizadores
  start <<- panel[c("tr","ts","a","n")]
  ## diagrama de dispersão
  plot(umid~log10(tens), data=da,
       ylab=expression(Conteúdo~de~água~(g~g^{-1})),
       xlab=expression(log[10]~Tensão~matricial~(kPa)),
       main=expression(f(x)==tr+(ts-tr)/(1+(a*x)^n)^{1-1/n}))
  ## sobrepõe a curva controlada pelos deslizadores
  with(start, curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n),
                    add=TRUE, col=2, lty=2))
  panel
}

par(mar=c(4.1,4.2,3.1,1))
# cria objetos vazios que serão preenchidos durante processo
da <- c(); start <- list(); aju <- list(); i <- c()
# abre interface gráfica
panel <- rp.control()
# controla parâmetros do modelo
rp.slider(panel, ts, 0.4, 0.8, initval=0.6, showvalue=TRUE, action=vg)
rp.slider(panel, tr, 0, 0.5, initval=0.3, showvalue=TRUE, action=vg)
rp.slider(panel, a, 0.1, 5, initval=1.3, showvalue=TRUE, action=vg)
rp.slider(panel, n, 1, 5, initval=1.6, showvalue=TRUE, action=vg)
# seleciona o conjunto de dados
rp.listbox(panel, caso, vals=levels(cras$caso),
           title="Subconjunto", action=vg)
# cria botão "Ajustar"
rp.button(panel, action=nlsajust, title="Ajustar")

sapply(aju, coef)    # estimativas dos coeficientes
lapply(aju, summary) # summary dos modelos ajustados

#-----------------------------------------------------------------------------
# dados de peso de capulhos em função da desfolha do algodão

cap <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt",
                  header=TRUE, sep="\t", encoding="latin1")
str(cap)
cap$desf <- cap$desf/100
cap <- subset(cap, select=c(estag, desf, pcapu))
cap$estag <- factor(cap$estag, labels=c("vegetativo","botão floral",
                                 "florescimento","maçã","capulho"))
str(cap)

xyplot(pcapu~desf|estag, data=cap, layout=c(5,1),
       xlab="Nível de desfolha artificial", ylab="Peso de capulhos")

#-----------------------------------------------------------------------------
# o exemplo a seguir mostra como usar com gráficos do pacote lattice

# função adapatada para outro modelo de regressão não linear
nlsajust <- function(panel){
  n0 <- try(nls(pcapu~f0-f1*desf^exp(C), data=da, start=start))
  if(class(n0)=="try-error"){
    trellis.focus("panel", nivel, 1, highlight=FALSE)
    grid.text(x=0.5, y=0.5, label="Não convergiu!\nAproxime mais.",
              gp=gpar(col="red"))
    trellis.unfocus()
  }
  trellis.focus("panel", nivel, 1, highlight=FALSE)
  with(as.list(coef(n0)), panel.curve(f0-f1*x^exp(C), add=TRUE, col=2))
  trellis.unfocus()
  print(coef(n0))
  aju[[i]] <<- n0
  panel
}

# função adapatada para outro modelo de regressão não linear
ptn <- function(panel){
  nivel <<- as.numeric(panel$nivel)
  i <<- levels(cap$estag)[nivel]
  da <<- subset(cap, estag==i)
  start <<- panel[c("f0","f1","C")]
  print({
    xyplot(pcapu~desf|estag, data=cap, layout=c(5,1), col=1,
           xlab="Nível de desfolha artificial", ylab="Peso de capulhos",
           main=expression(f(x)==f[0]-f[1]*x^exp(C)),
           strip=strip.custom(bg="gray90"))
  })
    trellis.focus("panel", nivel, 1, highlight=FALSE)
    with(start, panel.curve(f0-f1*x^exp(C), add=TRUE, col=2, lty=2))
    trellis.unfocus()
  panel
}

da <- c(); start <- list(); aju <- list(); i <- c(); nivel <- c()
panel <- rp.control()
rp.slider(panel, f0, 10, 50, initval=30, showvalue=TRUE, action=ptn)
rp.slider(panel, f1, 0, 30, initval=10, showvalue=TRUE, action=ptn)
rp.slider(panel, C, -5, 5, initval=0, showvalue=TRUE, action=ptn)
rp.listbox(panel, nivel, vals=1:nlevels(cap$estag),
           title="Painel", action=ptn)
rp.button(panel, action=nlsajust, title="Ajustar")

sapply(aju, coef)
lapply(aju, summary)

#-----------------------------------------------------------------------------

Ajuste do modelo potência à dados de redução na produção de algodão devido à desfolha com interface gráfica desenvolvida com as ferramentas do pacote rpanel.

Ajuste simultâneo de regressão/anova para várias respostas

Valores observados e ajustados acompanhados do intervalo de 95% em função da dose de potássio e nível de água para as variáveis resposta estudadas.

Na realização de experimento é comum a observação de várias variáveis respostas em uma mesma unidade amostral. Em termos estatísticos, o que observamos é um vetor aleatório. Por exemplo, em experimento de com soja para fins de ensaio de competição de cultivares o pesquisador pode observar: altura da planta, altura da primeira vagem, massa seca da planta, número de vagens por planta, teores foliares de N, P, K, Mg, Ca, peso de 100 grãos, proporção de sementes viáveis, dias para fechar o ciclo, todas elas observadas em cada unidade experimental. Nesse caso então, observamos um vetor aleatório de 12 variáveis aleatórias. O ideal seria empregarmos um modelo multivariado, porém como as respostas podem ser contínuas não limitadas, contagens, contínuas limitadas, contínuas com censura, proporções, etc, ainda não dispomos de um modelo de distribuição multivariado que contemple essa situação. O mais comum é fazer a análise univariada e como todas essas variáveis estão sob o efeito das mesmas fontes de variação, a análise estatística vai empregar o mesmo modelo (pelo menos no tocante a parte determinística). Sendo assim, tem como fazer o estudo de todas as variáveis reposta de uma vez só? Sim. Para saber como continue lendo.

Existem diversas formas de aplicar tarefas em grupos ou colunas. A família apply do stats oferece recurso para um número bem grande de situações. A família ply do plyr acrescenta mais opções. Além destes temos o doBy com mais algumas funcionalidades.

Bem, vamos ao que interessa. No CMR eu faço ajustes de regressão com polinômios para dados de um experimento com a cultura da soja (vasos em casa de vegetação). Foram estudados o nível de água (umidade do solo em relação à capacidade de campo) e potássio aplicado ao solo. O objetivo desse experimento foi verificar se o fornecimento de potássio pode minimizar impactos do deficit hídrico. Para acessar o artigo completo clique aqui: Umidade do solo e doses de potássio na cultura da soja. Não deixe de conferir a lista de links abaixo. Em próximos artigos do Rídiculas vamos aprofundar o estudo nesses dados comparando curvas. Até lá.

  • A brief introduction to “apply” in R;
  • plyr – The split-apply-combine strategy for R;
  • plyr: Tools for splitting, applying and combining data;
  • R Grouping functions: sapply vs. lapply vs. apply. vs. tapply vs. by vs. aggregate vs.;
  • Say it in R with “by”, “apply” and friends;
  • The doBy package;
  • doBy: doBy – Groupwise summary statistics, general linear contrasts, population means (least-squares-means), and other utilities;
  • #-----------------------------------------------------------------------------
    # importa dados
    soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                       header=TRUE, sep="\t", dec=",")
    str(soja)
    soja <- transform(soja, a=agua, A=factor(agua),
                      k=potassio, K=factor(potassio))
    resp <- 4:7 # índice das variáveis resposta
    str(soja)
    
    #-----------------------------------------------------------------------------
    # ver
    
    soja2 <- cbind(stack(soja[,resp]), soja[,-resp])
    
    require(lattice)
    
    xyplot(values~k|ind, groups=A, data=soja2,
           scales="free", jitter.x=TRUE, type=c("p","a"))
    
    #-----------------------------------------------------------------------------
    # ajuste do modelo saturado de efeitos categóricos (modelo de médias)
    
    m0 <- apply(soja[,resp], 2,
                function(r){ aov(r~bloco+A*K, soja) })
    lapply(m0, anova)
    
    #-----------------------------------------------------------------------------
    # ajuste de polinômio para potassio por nível de água
    
    m1 <- apply(soja[,resp], 2,
                function(r){ aov(r~bloco+A/poly(k, 4), soja) })
    lapply(m1, anova)
    names(coef(m1[[1]]))
    
    inter <- m1[[1]]$assign==3        # indices das interações
    eff <- names(coef(m1[[1]]))[inter]
    togrep <- outer(c("A37.5","A50","A62.5"), c("1$","2$","(3|4)$"),
                    paste, sep=":.*") # nomes de busca
    or <- order(togrep)
    tolist <- outer(c("A37.5","A50.0","A62.5"), c("lin","qua","lof"),
                    paste, sep="|")   # nomes de rótulo
    list.desd <- sapply(togrep, grep, eff)
    names(list.desd) <- tolist
    list.desd <- list.desd[or]
    list.desd                         # lista com os índices dos desdobramentos
    
    # anovas com desdobramento de efeito de potássio em cada nível de água
    lapply(m1, summary, split=list("A:poly(k, 4)"=list.desd))
    
    #-----------------------------------------------------------------------------
    # com base nas anovas acima selecionamos os modelos
    
    form <- list( rengrao~bloco+A/poly(k, 2, raw=TRUE),
                 pesograo~bloco+A/poly(k, 2, raw=TRUE),
                    kgrao~bloco+A/poly(k, 2, raw=TRUE),
                    pgrao~bloco+A/poly(k, 3, raw=TRUE))
    names(form) <- names(soja)[resp]
    
    m2 <- lapply(form, aov, data=soja, contrasts=list(bloco=contr.sum))
    lapply(m2, coef)
    
    #-----------------------------------------------------------------------------
    # predição a partir do modelo
    
    str(soja)
    pred <- expand.grid(bloco="I", A=levels(soja$A), k=seq(0,180,5))
    
    #aux <- lapply(m2, predict, newdata=pred, interval="confidence")
    aux <- lapply(m2, function(r){ # remove o efeito do bloco na predição
      predict(r, newdata=pred, interval="confidence")-coef(r)["bloco1"]})
    aux <- lapply(aux, as.data.frame)
    str(aux)
    
    require(plyr)
    pred <- cbind(ldply(aux), pred)
    pred$.id <- factor(pred$.id, levels=levels(soja2$ind))
    
    #-----------------------------------------------------------------------------
    # gráfico final
    
    fl <- c("K nos grãos (mg)", "Peso 100 grãos (g)",
            "P nos grãos (mg)", "Rendimento de grãos (g)")
    cbind(levels(pred$.id), levels(soja2$ind), fl)
    
    require(latticeExtra)
    display.brewer.all()
    mycol <- brewer.pal(6, "Set1")
    
    panel.ciH <- function(x, y, ly, uy, subscripts, ...){
      y <- as.numeric(y)
      x <- as.numeric(x)
      or <- order(x)
      ly <- as.numeric(ly[subscripts])
      uy <- as.numeric(uy[subscripts])
      panel.polygon(c(x[or], x[rev(or)]), c(ly[or], uy[rev(or)]),
                    col=1, alpha=0.05, border=NA)
      panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
      panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
      panel.xyplot(x, y, subscripts=subscripts, ...)
    }
    
    #png(file="f031.png", width=500, height=500)
    trellis.par.set(superpose.symbol=list(col=mycol),
                    superpose.line=list(col=mycol),
                    strip.background=list(col="gray90"))
    xyplot(values~k|ind, groups=A, data=soja2,
           scales="free", jitter.x=TRUE,
           xlab="Potássio no solo", ylab=NULL,
           strip=strip.custom(factor.levels=fl),
           auto.key=list(title="Níveis de água (%)", cex.title=1.2,
             columns=3, points=FALSE, lines=TRUE))+
      as.layer(xyplot(fit~k|.id, groups=A, pred, type="l", lwd=1.5,
                      ly=pred$lwr, uy=pred$upr,
                      panel.groups=panel.ciH,
                      panel=panel.superpose))
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # pode-se tentar comparar as curvas ajustando modelos apropriados sob H0
    # exemplo, será que as curvas vermelha e verde para kgrao e pesograo não
    # podem ser representados pela mesma curva?
    # as curvas azul e ver para pgrao diferem de uma reta horizontal?
    
    #-----------------------------------------------------------------------------
    
    Categorias:regressão Tags:, , ,

    Reparametrização do modelo quadrático para ponto crítico

    Perfil de log-verossimilhaça representada pelo módulo da raíz da função deviance para cada um dos parâmetros do modelo quadrático reparametrizado. Linhas tracejadas indicam os limites dos intervalos de confiança de 99, 95, 90, 80 e 50%. Padrões de perfil diferentes de um V indicam intervalos assimétricos.

    O emprego de funções polinomiais é muito frequente em ciências agrárias. O polinômio de segundo grau ou modelo quadrático é a escolha mais natural e imediata de modelo para quando a variável dependente apresentam uma tendência/relação não linear com a variável independente. Não raramente, o emprego desse modelo motiva também a estimação do ponto crítico (mínimo ou máximo). Para isso, o que se faz é estimar os parâmetros do modelo e via regras do cálculo obter a estimativa pontual do ponto crítico. O problema é que normalmente se necessita e não se faz, por desconhecimento de métodos, associação de incerteza à essa estimativa. É sobre isso que vou tratar agora: como associar incerteza ao ponto crítico de um modelo quadrático.

    O modelo quadrático tem essa equação

    m(x) = a+b\cdot x + c\cdot x^2,

    em que a é o intercepto, ou seja, E(Y|x=0) = ab é valor da derivada do modelo na origem, ou seja, d(m_0(x))/dx = b, e c mede a curvatura. Derivando o modelo, igualando a zero e resolvendo, obtemos que o ponto crítico é

    x_m = -b/(2\cdot c),

    portanto, x_m é uma função das estimativas dos parâmetros, e como tal, deve apresentar incerteza associada.

    Uma forma de associar incerteza é empregar o método delta para encontrar erros padrões. No R, o mesmo está disponível em msm::deltamethod(). Outra forma é usar um simulação boostratp e aplicar nas amostras bootstrap a relação acima. E o terceiro, e que eu mais gosto pois tem uma série de vantagens, é usar o perfil de log-verossimilhança.

    Para usar o perfil de verossimilhança eu tive que escrever o modelo de forma que o ponto crítico fosse um parâmetro dele. Reparametrizar um modelo nem sempre é uma tarefa fácil. Nesse caso o procedimento escrito é curto e simples de entender. Pórem eu demorei cerca de uma hora com tentativa-erro até usar a seguinte abordagem: imagine uma função que tenha ponto de máximo. Podemos aproximar essa função ao redor de um valor x_0 por série de Taylor. Para uma aproximação de segundo grau temos

    f(x) \approx f(x_0)+f'(x_0)(x-x_0)+0.5\cdot f''(x_0)(x-x_0)^2.

    Agora imagine que x_0 é um ponto crítico da função. O termo linear da aproximação é 0 no ponto crítico por definição. Se a função que estamos aproximando for perfeitamente quadrática, esse procedimento nos leva a reparametrização que desejamos. Assim, o modelo quadrático reparametrizado é

    n(x) = y_m+c\cdot (x-x_m)^2.

    em que x_m é o nosso parâmetro de interesse, o ponto crítico, y_m é o valor da função no ponto crítico, ou seja, E(Y|x=x_m) = y_m, c é a métade da segunda derivada da função no ponto crítico, ou seja, c=0.5\cdot f''(x_m), é uma curvatura. Então nos reparametrizamos o model quadrático para ter o ponto crítico como parâmetro. Nessa reparametrização obtvemos 3 parâmetros com interpretação simples. Em termos de interpretação, se c=0 não temos curvatura, portanto, não temos ponto crítico e devemos ajustar um modelo linear da reta. Se temos curvatura podemos estimar o ponto crítico e o valor da função no ponto crítico. Então imagine que você estuda doses de nutrientes para uma cultura, não é fundamental saber qual a dose que confere o máximo e qual a produção esperada nessa dose? Melhor que isso é poder associar incerteza à essas quantidades.

    O modelo reparametrizado, diferente do original, é um modelo não linear. Por isso, para estimar os parâmetros no R usaremos a função nls() e não a lm(). (Visite os posts passados para ver mais coisas relacionadas à modelos não lineares.) Para obter os intervalos de perfil de log-verossimilhança usaremos os métodos confint() e profile(). Consulte a documentação dessas funções para saber o que elas de fato fazem.

    No CMR abaixo eu ilustro a aplicação desses métodos. É recomendando ao leitor que compare os intervalos de confiança, que repita os procedimento simulando os dados com outros parâmetros, tamanhos e espaços, e se for o caso, que use dados reais. Futuros artigos do Rídiculas vão incluir comparação de modelos quadráticos (teste de hipótese) e o modelo quadrático para duas regressoras (bidimensional). Até a próxima rídicula.

    #-----------------------------------------------------------------------------
    # dados simulados
    
    set.seed(23456)
    x <- 1:12
    y <- 3+0.17*x-0.01*x^2+rnorm(x,0,0.05)
    plot(x, y)                         # dados observados
    curve(3+0.17*x-0.01*x^2, add=TRUE) # curva que gerou os dados
    
    #-----------------------------------------------------------------------------
    # ajuste do y = a+b*x+c*x², é um modelo linear
    
    m0 <- lm(y~x+I(x^2))
    summary(m0)
    
    #-----------------------------------------------------------------------------
    # método delta
    
    require(msm)
    coef <- coef(m0)[2:3]; vcov <- vcov(m0)[2:3,2:3]
    xm <- -0.5*coef[1]/coef[2]; xm                         # estimativa
    sd.xm <- deltamethod(~(-0.5*x1/x2), coef, vcov); sd.xm # erro padrão
    xm+c(-1,1)*qt(0.975, df=df.residual(m0))*sd.xm         # IC 95%
    
    #-----------------------------------------------------------------------------
    # simulação (bootstrap paramétrico)
    
    require(MASS)
    aa <- mvrnorm(1000, mu=coef, Sigma=vcov)
    xm.boot <- apply(aa, 1, function(i) -0.5*i[1]/i[2])
    mean(xm.boot)                      # média bootstrap
    sd(xm.boot)                        # desvio-padrão bootstrap
    quantile(xm.boot, c(0.025, 0.975)) # IC boostrap
    
    #-----------------------------------------------------------------------------
    # ajuste do y = ym+c*(x-xm)², é um modelo não linear
    
    plot(x, y)
    abline(h=3.7, v=8.5, lty=3) # aproximação visual do chute
    start <- list(ym=3.7, xm=8.5, c=coef(m0)[3])
    max(fitted(m0))             # boa aproximação de chute
    x[which.max(fitted(m0))]    # boa aproximação de chute
    n0 <- nls(y~ym+c*(x-xm)^2, start=start)
    summary(n0) # estimativa do xm que é parâmetro do modelo, estimado diretamente
    
    #-----------------------------------------------------------------------------
    # gráfico dos observados com a curva ajustada
    
    plot(x, y)
    with(as.list(coef(n0)),{
      curve(ym+c*(x-xm)^2, add=TRUE)
      abline(h=ym, v=xm, lty=3)
    })
    
    #-----------------------------------------------------------------------------
    # veja a igualdade dos modelos
    
    cbind(coef(m0),
          with(as.list(coef(n0)), c(ym+c*xm^2, -2*c*xm, c)))
    deviance(m0)
    deviance(n0)
    plot(fitted(m0), fitted(n0)); abline(0,1)
    
    #-----------------------------------------------------------------------------
    # intervalos de confiança
    
    confint.default(n0) # IC assintóticos, são simétricos por construção
    confint(n0)         # IC baseado em perfil de log-verossimilhança
    
    #-----------------------------------------------------------------------------
    # gráficos de perfil de verossimilhança
    
    #png(file="f030.png", width=500, height=175)
    par(mfrow=c(1,3))
    plot(profile(n0))
    layout(1)
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # representação do ajuste
    
    pred <- data.frame(x=seq(0,13,l=30))
    pred <- cbind(as.data.frame(predict(m0, newdata=pred, interval="confidence")),
                  pred)
    str(pred)
    ic <- confint(n0); ic
    
    plot(x, y, type="n")              # janela vazia
    with(pred,                        # bandas de confiança
         polygon(c(x,rev(x)), c(lwr,rev(upr)), col="gray90", border=NA))
    points(x, y)                      # valores observados
    with(pred, lines(x, fit))         # valores preditor
    abline(v=ic[2,], h=ic[1,], lty=2) # IC para xm e ym
    
    #-----------------------------------------------------------------------------
    

    Ilustração do método Newton-Raphson em regressão não linear

    Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação de parâmetros em um modelo de regressão não linear.

    Em regressão não linear, as estimativas dos parâmetros são obtidas por emprego de algum método numérico. Um dos métodos empregados é o de Newton-Raphson. Tal método é usado para encontrar raízes de uma função. Portanto, não é um método de maximação/mínimização de funções, mas usa o fato equivalente de que a derivada no ponto crítico de uma função é zero.

    Documentação sobre o método Newton-Raphson existe em abundância. O que motivou fazer esse post foi não ter encontrado nenhum que exemplificasse aplicação para estimação de parâmetros em um modelo de regressão não linear. O estimador de mínimos quadrados para um modelo de regressão é aquele que torna mínima a função

    S(\theta) = \sum_{i=1}^n (y_i-f(x_i, \theta))^2.

    O método Newton-Raphson vai procurar as raízes da derivada dessa função com relação à cada parâmetro. Então, se \theta é um vetor de p parâmetros, a derivada de S(\theta) será uma função p-dimensional,

    S_1(\theta) = \displaystyle\frac{\partial S(\theta)}{\partial\theta}=\begin{bmatrix}  \displaystyle\frac{\partial S(\theta)}{\partial\theta_1}\\  \vdots\\  \displaystyle\frac{\partial S(\theta)}{\partial\theta_p}  \end{bmatrix}.

    O método se baseia na aproximação de primeira ordem da função S_1(\theta) de forma que um valor da função para um valor \theta^{i+1} na redondeza de um valor \theta^i é

    S_1(\theta^{i+1}) \approx S_1(\theta^i)+ (\theta^{i+1}-\theta^i)\displaystyle\frac{\partial S_1(\theta)}{\partial\theta}|_{\theta=\theta^i}.

    O valor de S_1(\theta^{i+1}) é zero uma vez que \theta^{i+1} é a raíz da função considerando a expansão de primeira ordem, assim, reorganizando a expressão, temos que a atualização dos valores é data por

    \theta^{i+1} = \theta^i- \displaystyle\frac{S_1(\theta^i)}{S_2(\theta^i)},

    sendo que

    S_2(\theta) = \frac{\partial S_1(\theta)}{\partial\theta}=\begin{bmatrix}  \frac{\partial S_1(\theta)}{\partial\theta_1\partial\theta_1} & \ddots\\  \ddots& \frac{\partial S_1(\theta)}{\partial\theta_p\partial\theta_p}  \end{bmatrix}

    é a matriz de derivadas que tem dimensão p\times p.

    Esse procedimento numérico foi implementado no script abaixo para estimação dos parâmetros do modelo Michaelis-Menten e o modelo Exponencial reparametrizado que desenvolvi na minha Dissertação de Mestrado. No início eu monto as funções separadas e no final eu faço uma função que recebe os ingredientes mínimos e faz a estimação dos parâmetros. No final do script tem um procedimento gráfico para seleção dos chutes iniciais e acompanhar a trajetória durante as iterações. Até a próxima ridícula.

    #-----------------------------------------------------------------------------
    # método Newton-Raphson
    # objetivo é mínimizar a soma de quadrados dos resíduos mas iremos
    # encontrar as raízes da derivada dessa função pelo método de Newton-Raphson
    
    #-----------------------------------------------------------------------------
    # parâmetros, modelo e dados simulados deles
    
    theta <- c(A=10, B=2) # parâmetros
    x <- 1:12
    f <- function(theta, x){
      ## função que retorna os valore preditos para theta e x fornecidos
      ## theta: vetor de parâmetros A e B do modelo Michaelis-Menten
      ## x: vetor de valores da variável independente
      haty <- theta[1]*x/(theta[2]+x)
      return(haty)
    }
    set.seed(222); e <- rnorm(x,0,0.1)
    y <- f(theta=theta, x=x)+e
    plot(y~x)
    curve(f(theta, x), add=TRUE, col=2)
    
    #-----------------------------------------------------------------------------
    # função objetivo, derivada primeira e segunda com relação à theta
    
    S0 <- function(theta, f, y, x){
      ## função que retorna a soma de quadrados dos resíduos para
      ## theta: vetor de parâmetros do modelo f
      ## f: função do modelo de regressão não linear que depende de theta e x
      ## y e x: vetores dos valores das variáveis dependente e independente
      sqr <- sum((y-f(theta, x))^2)
      return(sqr)
    }
    
    S0(theta=theta, f=f, y=y, x=x)
    
    # derivada primeira, dimensão 1 x p
    D(expression((y-A*x/(B+x))^2), "A") # derivada de S0 com relação a A
    D(expression((y-A*x/(B+x))^2), "B") # derivada de S0 com relação a B
    
    # derivada segunda, dimensão p x p
    D(D(expression((y-A*x/(B+x))^2), "A"), "A") # derivada de S0 com relação a A e A
    D(D(expression((y-A*x/(B+x))^2), "B"), "B") # derivada de S0 com relação a B e B
    D(D(expression((y-A*x/(B+x))^2), "A"), "B") # derivada de S0 com relação a A e B
    D(D(expression((y-A*x/(B+x))^2), "B"), "A") # derivada de S0 com relação a B e A
    
    #-----------------------------------------------------------------------------
    # para ter essas funções vamos usar a deriv3()
    
    ff <- deriv3(~(y-A*x/(B+x))^2,
                 c("A","B"),
                 function(x, y, A, B){ NULL })
    
    ff(x=x, y=y, theta[1], theta[2])
    
    #-----------------------------------------------------------------------------
    # mas as matrizes são 1 x p e p x p, então falta fazer a soma em n
    
    fff <- function(x, y, A, B){
      ## função que retorna gradiente e hessiano avaliado em theta
      ## y e x: vetores dos valores das variáveis dependente e independente
      ## A e B: valores dos parâmetros do modelo Michaelis-Menten
      m <- ff(x, y, A, B)
      g <- attr(m, "gradient"); G <- apply(g, 2, sum)
      h <- attr(m, "hessian"); h <- matrix(h, ncol=2*2)
      H <- apply(h, 2, sum); H <- matrix(H, 2, 2)
      return(list(G=G, H=H))
    }
    
    fff(x=x, y=y, theta[1], theta[2])
    
    #-----------------------------------------------------------------------------
    # agora é só fazer o Newton-Raphson
    
    th0 <- c(9,2); tol <- 1e-10; niter <- 100;
    dis <- 1; i <- 0
    while(dis>tol & i<=niter){
      GH <- fff(x=x, y=y, th0[1], th0[2])
      G <- GH[[1]]
      iH <- solve(GH[[2]])
      th1 <- th0-iH%*%G
      dis <- sum((th0-th1)^2)
      i <- i+1
      print(c(i, dis, th1))
      th0 <- th1
    }
    
    n0 <- nls(y~A*x/(B+x), start=c(A=9, B=2))
    coef(n0)
    
    iH/vcov(n0) # hessiano é proporcional à matriz de covariância dos parâmetros
    
    #-----------------------------------------------------------------------------
    # vamos incluir todos os elementos em uma única função!
    
    newton <- function(formula, data, start, tol=1e-10, niter=100, trace=TRUE){
      ## formula: expressão da resposta como função das covariáveis e parâmetros
      ## data: data.frame que contém dos dados
      ## start: vetor ou lista nomeados de valores iniciais para os parâmetros
      ## tol: norma mínima entre valores consecutivos dos parâmetros
      ## niter: número máximo de iterações do algortimo
      ## trace: se TRUE imprime o histórico de iterações
      form <- paste("~(",paste(formula[2:3], collapse="-"),")^2", sep="")
      vars <- all.vars(formula)
      nvars <- length(vars)
      pars <- 1:(nvars+1)
      l <- as.list(pars)
      names(l)[1:nvars] <- vars
      f <- as.function(l)
      body(f) <- NULL
      ff <- deriv3(as.formula(form), names(start), f)
      th0 <- unlist(start)
      i <- 0; dis <- 100; iter <<- matrix(th0, ncol=length(start))
      while(dis>tol & i<=niter){
        m <- do.call(ff, c(as.list(data), as.list(th0)))
        g <- attr(m, "gradient"); G <- apply(g, 2, sum)
        h <- attr(m, "hessian"); h <- matrix(h, ncol=length(start)^2)
        H <- apply(h, 2, sum); H <- matrix(H, nrow=length(start))
        iH <- solve(H)
        th1 <- th0-iH%*%G
        dis <- sum((th0-th1)^2)
        i <- i+1
        iter <<- rbind(iter, t(th1))
        if(trace==TRUE) print(c(i, dis, th1))
        th0 <- th1
      }
      th0 <- as.vector(th0)
      names(th0) <- names(start)
      ## retorna o último valor iterado e a uma matriz com o histórico
      return(list(theta=th0, iter=iter))
    }
    
    #-----------------------------------------------------------------------------
    # vamos estudar vetor de chutes presentes nessa lista
    
    start.l <- list(c(A=10,B=2.5), c(A=10,B=0.1), c(A=14,B=2.5),
                    c(A=15.13, B=1.7), c(A=18.4,B=0.13))
    
    data <- data.frame(x=x, y=y)
    formula <- y~A*x/(B+x)
    start <- start.l[[4]]
    nr <- newton(formula, data, start)
    nr$theta
    
    sqe <- function(A, B, y, x){ hy <- (A*x)/(B+x); sum((y-hy)^2) }
    SQE <- Vectorize(sqe, c("A", "B"))
    
    A.grid <- seq(-3,15,l=100)
    B.grid <- seq(-3,15,l=100)
    sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
    #png("f023.png", w=500, h=500)
    contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
            xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
    for(i in seq(nrow(iter)-1)){
      arrows(iter[i,1], iter[i,2],
             iter[i+1,1], iter[i+1,2],
             col=2, length=0.1)
      Sys.sleep(0.1)
    }
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # com o mouse selecione regiões da superfície para serem valores iniciais
    
    loc <- locator()
    loc <- as.data.frame(loc)
    names(loc) <- c("A","B")
    start.l <- split(loc, f=1:nrow(loc))
    
    #-----------------------------------------------------------------------------
    # ver como fica com o modelo Exponencial
    
    n0 <- nls(y~A*(1-exp(-log(2)*x/B)), start=c(A=9, B=2))
    coef(n0)
    
    data <- data.frame(x=x, y=y)
    formula <- y~A*(1-exp(-log(2)*x/B))
    start <- list(A=5, B=0.5)
    nr <- newton(formula, data, start)
    nr$theta
    
    sqe <- function(A, B, y, x){ hy <- A*(1-exp(-log(2)*x/B)); sum((y-hy)^2) }
    SQE <- Vectorize(sqe, c("A", "B"))
    
    A.grid <- seq(-10,20,l=100)
    B.grid <- seq(-10,20,l=100)
    sqe.surf <- outer(A.grid, B.grid, SQE, data$y, data$x)
    contour(A.grid, B.grid, sqe.surf, levels=(1:45)^2, asp=1,
            xlab=expression(theta[0]), ylab=expression(theta[1]), col="gray70")
    for(i in seq(nrow(iter)-1)){
      arrows(iter[i,1], iter[i,2],
             iter[i+1,1], iter[i+1,2],
             col=2, length=0.1)
      Sys.sleep(0.1)
    }
    #-----------------------------------------------------------------------------
    
    

    Como fazer regressão polinomial com diferentes graus

    Ajuste de modelos de regressão polinomial com diferentes graus para o polinômio. Para as cultivares da esquerda o polinômio é quadrático, para a da direita o polinômio é cúbico.

    No post Predição em modelos de regressão com blocos foram apresentados procedimentos para estimação e predição de parâmetros em modelos de regressão com blocos. Nesse post aqui eu vou estender o exemplo permitindo o ajuste com diferentes graus para o polinômio que descreve a relação da resposta com a preditora contínua. Nos códigos abaixo vou apresentar aspectos dos polinômios ortogonais, do cálculo do R² e de como fazer o ajuste com diferentes graus de polinômio em um mesmo modelo, ou em outras palavras, em uma mesma lm(). No final, como de costume, confecciono o gráfico resultado do modelo ajustado com as bandas de confiança. Mais uma vez deixei o código comentado passo a passo para orientar o leitor. Até a próxima ridícula!

    #-----------------------------------------------------------------------------
    # dados de índice agronômico de milho de um experimento
    
    da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                     header=TRUE, sep="\t")
    str(da)
    
    #-----------------------------------------------------------------------------
    # objetivo: apresentar procedimentos para ajuste de polinômios e como ajustar
    # polinômios de diferentes graus em uma mesma análise
    
    #-----------------------------------------------------------------------------
    # vamos apresentar essas funções para o caso do ajuste de uma função e fazer
    # de conta que não temos blocos apenas por facilidade de exposição
    
    levels(da$cultivar)
    br <- subset(da, cultivar=="BR-300")
    str(br)
    
    #-----------------------------------------------------------------------------
    # ajustar diminuindo o nível do polinômio até que não tenhamos falta de ajuste
    # fazer isso de diferentes formas:
    # 1) fórmula explicita e anova seqüencial
    # 2) polinômios ortogonais e teste t
    
    # 1) fórmula explicita e anova seqüencial
    m1.5 <- lm(indice~dose+I(dose^2)+I(dose^3)+I(dose^4)+I(dose^5), data=br)
    anova(m1.5) # testes F apontam que os termos > 2 grau não contribuem
    m1.2 <- lm(indice~dose+I(dose^2), data=br)
    anova(m1.2, m1.5) # o abandono desses termos não resulta em falta de ajuste
    anova(m1.2)
    
    # 2) polinômios ortogonais e teste t
    m2.5 <- lm(indice~poly(dose, degree=5), data=br)
    anova(m2.5)   # anova não separa os termos
    summary(m2.5) # quando usamos polinômios orotoganis, o valor t^2 = F
    cbind("t^2"=summary(m2.5)$coeff[-1,3]^2, "F"=anova(m1.5)[-6,4]) # viu?
    
    m2.2 <- lm(indice~poly(dose, degree=2), data=br)
    summary(m2.2)
    
    #-----------------------------------------------------------------------------
    # antigamente experimentos agronômicos com fatores quantitativos eram
    # planejados para polinômios ortogonais. assim usava espaçamento regular entre
    # níveis. muitas pessoas acham que isso é regra mas era apenas algo imposto
    # pelas limitações de computação da época. com polinômios ortogonais a matriz
    # X'X é ortogonal, então facilmente invertível e obtém-se as estimativas dos
    # parâmetros. com essas estimativas e erros padrões aplicava-se o teste t (ao
    # invés do teste F) para testar a significância dos termos. os temos não
    # significativos eram abandonados mas a estimação não precisava ser refeita
    # por causa da ortogonalidade. observe os procedimentos abaixo.
    
    X <- model.matrix(m2.5) # matriz de delineamento ortogonal
    round(colSums(X), 4) # soma nas colunas dá zero
    round(t(X)%*%X, 4)   # ortogonal, portanto a inversa é a inversa da diagonal
    t(X)%*%br$indice     # as estimativas dos parâmetros, intercepto dividir por n
    
    cbind(X[,2:3], model.matrix(m2.2)[,2:3])[1:6,] # iguais colunas
    coef(m2.5); coef(m2.2) # as estimativas são iguais, não requer outro ajuste
    
    #-----------------------------------------------------------------------------
    # no entanto, esses coeficientes não tem interpretação pois a matriz não está
    # na escala das doses. aliás, polinômio não tem interpretação em nenhuma
    # escala. eles representam apenas uma aproximação local. interpretar a
    # magnitude das estimativas é perda de tempo. compare apenas o rank e a
    # direção (sinal). faça o teste de hipótese e interprete os valores preditos.
    # caso queira apresentar a equação ajustada, terá que obter as estimativas
    # sem polinômios ortogonais.
    
    #-----------------------------------------------------------------------------
    # obter as estimativas sem polinômios ortogonais com a função poly()
    
    m2.2 <- lm(indice~poly(dose, degree=2, raw=TRUE), data=br)
    summary(m2.2) # estimativas na escala do fator, são as mesmas do modelo m1.2
    
    #-----------------------------------------------------------------------------
    # muitos tutoriais/aplicativos obtém o R² do ajuste fazendo ajuste às médias.
    # *eu* não considero essa prática porque esse R² não reflete a razão
    # ruído/sinal dos dados para com o modelo, mas da média dos dados para com o
    # modelo. esse tipo de procedimento só é possível quando temos repetição de
    # níveis. veja os três exemplos simulados a seguir
    
    desv <- rnorm(12,0,1)      # desvios aleatórios
    a <- 0; b <- 1             # valores dos parâmetros
    x1 <- seq(1,12,1)          # sem repetição de nível
    x2 <- rep(seq(3,12,3),e=3) # com repetição de nível
    y1 <- a+x1*b+desv          # valores observados de y1
    y2 <- a+x2*b+desv          # valores observados de y2
    
    m1 <- lm(y1~x1); summary(m1) # R² correto
    m2 <- lm(y2~x2); summary(m2) # R² correto
    
    # R² calculado com relação ao desvio dos ajustados para as médias
    md <- tapply(y2, x2, mean)            # tiramos a média das repetições
    aj <- unique(round(fitted(m2),5))     # valores ajustados pelo modelo
    r2 <- 1-sum((md-aj)^2)/sum((md-mean(md))^2); r2 # R² como desvio das médias
    
    #-----------------------------------------------------------------------------
    # não sei dizer qual a justificativa para o uso desse procedimento para obter
    # o R². é um procedimento atraente porque o R² será *sempre* maior e talvez
    # essa seja a razão da sua adoção (que na *minha opinião* é errada). esse R²
    # não associa desvio das observações. então o valor desse R² não é a
    # porcentagem de explicação que o modelo dá para as observações, e sim para a
    # média delas! como assim se estou modelando as observações e não a média
    # amostral delas? veja esse exemplo mais abrupto
    
    #-----------------------------------------------------------------------------
    # criando um novo par x, y
    
    x3 <- x2+c(-0.001,0,0.001) # deslocamento mínimo faz perder as repetições
    y3 <- a+x3*b+desv
    
    cbind(y2,y3,x2,x3)
    
    # ambos pares x, y são praticamente a mesma coisa e o R² será muito próximo
    # para ambos. porém, para x2 pode-se usar o R² relativo a desvio da média.
    # se eu calcular esse R² ele vai diferir muito pois representa outra razão
    # ruído/sinal. *eu* não gosto disso.
    
    m3 <- lm(y3~x3)
    c(summary(m1)$r.sq, r2, summary(m2)$r.sq, summary(m3)$r.sq) # compara R²
    
    #-----------------------------------------------------------------------------
    # para mais críticas sobre R² procure sobre a análise do quarteto de Anscombe.
    
    #-----------------------------------------------------------------------------
    # ajustar modelos de regressão para os 3 níveis de cultivar e decompor a SQ.
    # aqui usaremos o fator de níveis ordenados que tem como padrão o contraste
    # de polinômios ortogonais. na poly() conseguimos especificar o grau e no
    # fator ordenado não conseguimos, mas as somas de quadrados podem ser
    # facilmente desdobradas se usarmos aov().
    
    m0 <- aov(indice~bloco+cultivar*ordered(dose), data=da)
    summary(m0)    # pode-se usar anova(m0) também
    summary.lm(m0) # quadro de estimativas dos parâmetros
    
    #-----------------------------------------------------------------------------
    # para obter o quadro de análise de variância com desdobramentos dos termos
    # do polinômio dentro dos níveis de cultivar precisamos mudar a fórmula do
    # modelo, é apenas uma mudança a nível de colunas mas o espaço do modelo ainda
    # é o mesmo. compare as colunas das matrizes obtidas com os comando abaixo.
    
    aux <- expand.grid(cultivar=gl(2,1), dose=1:3) # estrutura experimental
    model.matrix(~cultivar*ordered(dose), aux)     # cultivar cruzado com dose
    model.matrix(~cultivar/ordered(dose), aux)     # dose aninhado em cultivar
    
    #-----------------------------------------------------------------------------
    # ajustando modelo com fórmula corresponde ao desdobramento desejado
    # o procedimento abaixo é para criar a lista que especifica a partição nas
    # somas de quadrados, envolve algumas funções que vale a pena consultar o help
    
    da$D <- ordered(da$dose) # cria uma nova coluna apenas por conforto
    m1 <- aov(indice~bloco+cultivar/D, data=da) # modelo com a fórmula conveniente
    
    m1$assign # valores iguais a 3 representam os estimativas do 3 termo
    tm <- names(coef(m1))[m1$assign==3]                  # nomes dos termos
    tm <- gsub("^cultivar(.*):.*(.{2})$", "\\1-\\2", tm) # edita os nomes
    ord <- c(sapply(levels(da$cultivar), grep, x=tm, simplify=TRUE)) # ordena
    split <- as.list(ord)   # cria a lista com a posição dos termos
    names(split) <- tm[ord] # atribui nome aos elementos da lista
    
    #-----------------------------------------------------------------------------
    # quadro de análise de variância com somas que quadrados particionadas.
    # com ela escolhemos o grau do polinômio adequado para cada cultivar.
    
    summary(m1, split=list("cultivar:D"=split))
    
    #-----------------------------------------------------------------------------
    # uma das cultivares requer 3 grau, as outras ficam com o 2 grau. juntar os
    # termos não significativos do polinômio para criar o termo falta de ajuste.
    
    do.call(c, split) # ver os índices e fazer o agrupamento gerando nova split
    split <- list(Ag.L=1, Ag.Q=4,         Ag.lof=c(7,10,13),
                  BR.L=2, BR.Q=5,         BR.lof=c(8,11,14),
                  Pi.L=3, Pi.Q=6, Pi.C=9, Pi.lof=c(12,15))
    summary(m1, split=list("cultivar:D"=split)) # lof não significativos, ok!
    
    #-----------------------------------------------------------------------------
    # não iremos apresentar as estimativas do modelo m1 nem seus valores ajustados
    # porque esses se referem ao modelo completo. temos que ajustar o modelo
    # reduzido removendo os temos apontado pelo teste F. este modelo é o modelo
    # que apresentaremos as estimativas e faremos predição. note que o grau do
    # polinômio agora muda com o nível de cultivar. teremos que arrumar um
    # artifício ajustar um modelo considere isso. então vou criar uma indicadora
    # para multiplicar o termo cúbico e sumir com ele quando não for necessário.
    
    levels(da$cultivar)
    da$ind <- 0; da[da$cultivar=="Pioneer-B815",]$ind <- 1 # a indicadora do nível
    da$ind # o valor 1 é associado ao nível de cultivar que tem o termo cúbico
    
    #-----------------------------------------------------------------------------
    # o termo cúbico só vai existir para o nível Pionner por causa da indicadora.
    # vai aparecer um NA para os outros níveis
    
    m2 <- lm(indice~bloco+cultivar*(dose+I(dose^2)+I(ind*dose^3)), data=da)
    summary(m2)
    
    #-----------------------------------------------------------------------------
    # os NA implicam em algumas conseqüência como no funcionamento da predict().
    # para evitá-los podemos criar uma coluna para representar esse termo cúbico e
    # e declara uma fórmula de modelo que não retorne NA.
    
    da$Pi.C <- with(da, ind*dose^3) # é a coluna do cubo da dose para Pioneer
    m3 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+Pi.C, data=da,
             contrast=list(bloco=contr.sum))
    summary(m3)
    
    #-----------------------------------------------------------------------------
    # faz o gráfico final de predição pelo método matricial
    
    # data.frame da predição
    pred <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,10))
    pred$Pi.C <- with(pred, ifelse(cultivar=="Pioneer-B815", dose^3, 0))
    
    # matriz de predição sem colunas para os efeitos de blocos
    X <- model.matrix(~cultivar*(dose+I(dose^2))+Pi.C, data=pred)
    head(X)
    
    # vetor de estimativas sem os efeitos de blocos
    b <- coef(m3)[-grep("bloco", names(coef(m3)))]
    
    # vetor de erros padrões das estimativas
    se <- predict(m3, newdata=cbind(bloco=levels(da$bloco)[1], pred),
                  se.fit=TRUE)$se
    
    # quantil superior da distribuição t para um IC de 95%
    tc <- qt(0.975, df=df.residual(m3))
    
    # margem de erro
    me <- se*tc
    
    # data.frame com valores preditos e IC
    pred$fit <- X%*%b
    pred$lwr <- pred$fit-me; pred$upr <- pred$fit+me
    
    require(reshape)
    
    pred <- melt(pred[,-3], id=c("cultivar","dose"))
    names(pred)[ncol(pred)] <- "indice"
    da$variable <- "obs"
    pred <- rbind(da[,names(pred)], pred)
    
    #-----------------------------------------------------------------------------
    # gráfico com observados, preditos e IC pelo método matricial
    
    require(lattice)
    
    #png("f019.png", w=500, h=250)
    xyplot(indice~dose|cultivar, groups=variable, data=pred,
           distribute.type=TRUE, layout=c(3,1),
           type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2),
           strip=strip.custom(bg="gray90"),
           xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
           sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
           key=list(space="top", columns=3, type="o", divide=1,
             lines=list(pch=c(1,NA,NA), lty=c(0,1,2), col=c(1,2,3)),
             text=list(c("valores observados","valores preditos","IC 95%"))))
    #dev.off()
    
    #-----------------------------------------------------------------------------
    

    Predição em modelos de regressão com blocos

    Padrão de reposta de cultivares de milho à adubação com nitrogênio.

    Em experimentos onde não se dispõe de unidades experimentais uniformes (UI), usualmente se agrupa UI semelhantes em grupos que usualmente recebe o nome de bloco. Em cada bloco as unidades experimentais são uniformes, e entre classes são diferentes. Uniformidade aqui significa que entre essas UI não é possível atribuir uma causa/origem as pequenas diferenças existentes. São exemplos de blocos em experimentos agrícolas: nível topográfico, idade inicial, posição na casa de vegetação, operador do sistema, local de origem, horário de avaliação, dentre muitos outros que podem ser citados.

    O uso de um fator de blocagem permite que parâmetros e contrastes sejam estimados com maior precisão porque parte da diferença entre UI foi identificada pelo arranjo em blocos. No entanto, o uso de blocos sem a necessidade pode elevar a taxa de erro tipo I e o uso de blocos de forma incorreta pode reduzir o poder dos testes de hipótese. Ao iniciar uma pesquisa, consulte um especialista em planejamento de experimentos.

    Os blocos representam uma causa de variação presente mas que não é o foco da pesquisa. O pesquisador deve declarar o efeito de bloco no modelo estatístico da análise mas as inferências só serão aplicadas aos demais fatores do estudo. Em algum ponto, além de estimar e comparar parâmetros, há o interesse de fazer a predição da resposta. Nessa predição, normalmente, deseja-se que estejam presentes apenas os termos de interesse da pesquisa e remova-se o efeito dos blocos.

    Neste post vou apresentar procedimentos para obter valores preditos livre dos efeitos dos blocos. Farei uso da função predict() e mostrarei como fazer por operação de matrizes por meio de funções lineares. Deve-se fazer algumas adaptações nos procedimentos quando o experimento possui observações não disponíveis (parcelas perdidas). Algumas adaptações serão necessárias também para uso modelos lineares generalizados. Nesses procedimentos o efeito de bloco é fixo e não aleatório.

    Os dados são de experimento conduzido com parcelas dispostas em blocos. Os fatores de interesse são a cultivar de milho e a dose de nitrogênio. O objetivo da análise é obter o padrão de resposta de cada cultivar à adubação com nitrogênio. Os resultados serão apresentados por meio de gráfico. O código está comentado em casa sessão para orientar o leitor. Até a próxima ridícula.

    #-----------------------------------------------------------------------------
    # dados de índice agronômico de milho de um experimento
    
    da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                     header=TRUE, sep="\t")
    str(da)
    
    #-----------------------------------------------------------------------------
    # estrutura dos dados, experimento balanceado ortogonal
    # 4 níveis de bloco (identifica diferenças/variações de ambiente e manejo),
    # 3 níveis de cultivar de milho (fator de interesse)
    # 6 doses de nitrogênio igualmente espaçadas (fator de interesse)
    # 72 observações correspondentes ao produto cartesiano dos níveis dos fatores
    
    #-----------------------------------------------------------------------------
    # contraste polinomial, cria matriz de delineamento com colunas ortogonais.
    # para dois conjuntos de valores diferentes, o contr.poly é o mesmo
    # pois é estabelecido pelo rank e não pelos valores em si.
    # só é interpretável para seqüências regulares de valores.
    # isso é só demonstração, não iremos usar ortogonais nos modelos reduzidos.
    
    x <- ordered(seq(1,10,1))                     # seqüência regular
    y <- ordered(c(1,1.5,2.5,3,3.5,4,4.5,6,8,10)) # seqüência irregular
    X <- model.matrix(~x)    # matriz do modelo para níveis de x
    Y <- model.matrix(~y)    # matriz do modelo para níveis de y
    cbind(X[,2:3], Y[,2:3])  # observe que as matrizes são iguais
    round(t(X)%*%X)          # colunas ortogonais implica em matriz de covariância
                             # diagonal para as estimativas dos parâmetros
    
    #-----------------------------------------------------------------------------
    # objetivo: ajustar modelos de regressão para cada cultivar em função da dose,
    # modelos de regressão polinomial, no máximo de grau 3, fazer a predição livre
    # do efeito de bloco pois não há interesse em diferenças entre níveis.
    
    #-----------------------------------------------------------------------------
    # cria coluna de fator com níveis ordenados para dose
    
    unique(da$dose)           # doses únicas usadas
    da$ds <- ordered(da$dose) # fator ordenado
    str(da$ds)                # mostra a ordem dos níveis
    
    #-----------------------------------------------------------------------------
    # ajuste do modelo saturado, ou seja, grau máximo permitido para o polinômio
    # número de parâmetros estimados para dose é número de níveis - 1
    
    m0 <- lm(indice~bloco+cultivar*ds, data=da)
    par(mfrow=c(2,2)); plot(m0); layout(1) # checa as pressuposições
    
    #-----------------------------------------------------------------------------
    # ajuste do modelo de regressão linear, anova() verifica falta de ajuste (lof)
    
    m1 <- lm(indice~bloco+cultivar*dose, data=da)
    par(mfrow=c(2,2)); plot(m1); layout(1) # gráfico 1,1 aponta que requer termos
    anova(m0, m1)                          # compara modelos aninhados, testa lof
    
    #-----------------------------------------------------------------------------
    # ajusta o modelo quadrático para dar curvatura, 3 alternativas para isso
    # 1) fórmula explícita;
    # 2) polinômios ortogonais de grau k=degree;
    # 3) polinômios não ortogonais, assim, estimativas iguais em 1 e 3;
    # os valores ajustados/preditos não mudam com as alternativas.
    
    m2.1 <- lm(indice~bloco+cultivar*(dose+I(dose^2)), data=da)
    m2.2 <- lm(indice~bloco+cultivar*poly(dose, degree=2), data=da)
    m2.3 <- lm(indice~bloco+cultivar*poly(dose, degree=2, raw=TRUE), data=da)
    sapply(list(m2.1, m2.2, m2.3), coef)   # coeficientes lado a lado
    sapply(list(m2.1, m2.2, m2.3), fitted) # valores ajustados lado a lado
    
    m2 <- m2.1
    par(mfrow=c(2,2)); plot(m2); layout(1) # ok
    anova(m0, m2)                          # modelo quadrático não apresenta lof
    anova(m2)                              # teste de hipóteses com o modelo final
    
    #-----------------------------------------------------------------------------
    # como fazer predição
    # predição com fitted(): faz a predição total, considerando todos os termos do
    #                        modelo, tem dimensão dos dados, chama-se de valores
    #                        ajustados.
    # predição com predict(): faz a predição para quaisquer valores, mas todos os
    #                         termos estão presentes.
    # predição com produto de matrizes: faz o que você quiser.
    
    #-----------------------------------------------------------------------------
    # valores observados (o) e ajustados (f)
    
    of <- cbind(o=da$indice, f=fitted(m2)); of
    plot(of)
    
    #-----------------------------------------------------------------------------
    # usando a função predict() para prever usando apenas um nível de bloco.
    # vou usar o primeiro nível. como bloco ortogonal aos demais efeitos, os
    # contrastes entre quaisquer outros níveis dos fatores são os mesmos
    # independente do nível do bloco. fazer uma sequencia mais fina para a dose.
    
    bl <- levels(da$bloco)[1] # nível de bloco usando na predição
    ct <- levels(da$cultivar) # todos os níveis de cultivar serão usados
    ds <- seq(0,300,by=2)     # seqüência fina de valores para dose
    
    pred <- expand.grid(bloco=bl, cultivar=ct, dose=ds) # data.frame da predição
    aux <- predict(m2, new=pred, interval="confidence") # IC 95%
    str(aux) # matriz com colunas fit, lwr, upr
    
    #-----------------------------------------------------------------------------
    # tratamento dos dados para gráfico com bandas usando lattice::xyplot()
    
    pred1 <- cbind(pred, as.data.frame(aux))
    str(pred1) # contém apenas o nível I de bloco
    
    require(reshape)
    
    pred1 <- melt(pred1, id=c("cultivar","dose","bloco"))
    names(pred1)[ncol(pred1)] <- "indice"
    str(pred1)
    da$variable <- "obs"
    str(da)
    
    pred1 <- rbind(da[,c("cultivar","dose","bloco","variable","indice")], pred1)
    str(pred1) # junção dos valores observados e preditos
    
    #-----------------------------------------------------------------------------
    # gráfico com observados, preditos e IC
    
    require(lattice)
    
    xyplot(indice~dose|cultivar, groups=variable, data=pred1,
           distribute.type=TRUE, sub="predição para o bloco I",
           type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2))
    
    #-----------------------------------------------------------------------------
    # o gráfico aponta desvio sistemático dos observados para os preditos porque
    # nossa predição considerou só o bloco I. se mudarmos de bloco as curvas serão
    # movimentadas verticalmente pela mudança do intercepto para o de outro bloco.
    # para que o predito passe no meio dos dados é preciso fazer o efeito de bloco
    # se anular. se usarmos o contraste soma zero para blocos, podemos saber
    # quanto o bloco I difere de zero e sutrair esse valor do fit, lwr e upr
    
    #-----------------------------------------------------------------------------
    # usar contraste soma zero e diminuir dos preditos o efeito efeito do bloco I
    
    m3 <- lm(formula(m2), data=da, contrasts=list(bloco=contr.sum))
    cbind(coef(m2), coef(m3)) # os coeficientes de bloco mudaram
    
    blI <- coef(m3)["bloco1"]; blI # efeito do bloco I sob a restrição soma zero
    
    # remove o efeito do bloco I dos dados
    pred1[pred1$variable!="obs",]$indice <- pred1[pred1$variable!="obs",]$indice-blI
    
    #-----------------------------------------------------------------------------
    # faz o gráfico sem o efeito de bloco
    
    xyplot(indice~dose|cultivar, groups=variable, data=pred1,
           distribute.type=TRUE, sub="valores preditos sem efeito de bloco",
           type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2))
    
    #-----------------------------------------------------------------------------
    # predição pelo procedimento matricial. é um método mais geral. é valido para
    # modelos da classe glm() com algumas modificações.
    
    # matriz de predição sem colunas para os efeitos de blocos
    X <- model.matrix(~cultivar*(dose+I(dose^2)), data=pred)
    
    # vetor de estimativas sem os efeitos de blocos
    b <- coef(m3)[-grep("bloco", names(coef(m3)))]
    
    # vetor de erros padrões das estimativas
    se <- predict(m3, newdata=pred, se.fit=TRUE)$se
    
    # quantil superior da distribuição t para um IC de 95%
    tc <- qt(0.975, df=df.residual(m3))
    
    # margem de erro
    me <- se*tc
    
    # data.frame com valores preditos e IC
    pred2 <- pred[,-1]
    pred2$fit <- X%*%b
    pred2$lwr <- pred2$fit-me; pred2$upr <- pred2$fit+me
    pred2 <- melt(pred2, id=c("cultivar","dose"))
    names(pred2)[ncol(pred2)] <- "indice"
    
    pred2 <- rbind(da[,c("cultivar","dose","variable","indice")], pred2)
    str(pred2) # data.frame dos valores preditos pelo método matricial
    
    #-----------------------------------------------------------------------------
    # gráfico com observados, preditos e IC pelo método matricial
    
    #png("f018.png", w=500, h=250)
    xyplot(indice~dose|cultivar, groups=variable, data=pred2,
           distribute.type=TRUE, layout=c(3,1),
           type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2),
           strip=strip.custom(bg="gray90"),
           xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
           sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
           key=list(space="top", columns=3, type="o", divide=1,
             lines=list(pch=c(1,NA,NA), lty=c(0,1,2), col=c(1,2,3)),
             text=list(c("valores observados","valores preditos","IC 95%"))))
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # compara os métodos de predição
    
    round(pred1$indice-pred2$indice, 3) # pred1 e pred2 são iguais
    
    #-----------------------------------------------------------------------------
    # nesse procedimento consideramos os efeitos de bloco como fixos.
    # caso fossem de efeito aleatório poderiamos usar a função nlme::lme() e
    # declarar o efeito de blocos como aleatório. os valores preditos seriam
    # obtidos com predict(..., level=0). algumas manipulações seriam necessárias
    # para obter o intervalo de confiança.
    
    #-----------------------------------------------------------------------------
    

    Como adicionar legendas em gráficos da lattice

    Gráfico com bandas de confiança, legenda e assinatura de autor e versão.

    No post Como fazer legendas em gráficos eu mostrei procedimentos para gráficos do pacote graphics. Nesse post vou apresentar procedimento para fazer legendas em gráficos da lattice, desenvolvido por Deepayan Sarkar. Iremos empregar as opções auto.key=, key=, colorkey=, e as funções draw.key() e draw.colorkey(). O pacote grid, desenvolvido por Paul Murrell será necessário para posicionarmos as legendas na página gráfica. Nos exemplos temos legendas para pontos, linhas, pontos e linhas, barras, escala de cores para gráficos de três dimensões, gráfico com duas legendas, posicionamento da legenda com o mouse e gráfico com legenda para representar o ajuste de um modelo de regressão com valores observados, preditos e intervalo de confiança.

    #-----------------------------------------------------------------------------
    # pacotes
    
    require(lattice) # para fazer os gráficos
    require(grid)    # para posicionar as legendas
    
    #-----------------------------------------------------------------------------
    # dados de indice agronômico em função da cultivar de milho e dose de N
    
    da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                     header=TRUE, sep="\t")
    str(da)
    
    #-----------------------------------------------------------------------------
    # legenda para gráfico de dispersão
    
    # gráfico com legenda padrão
    xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE, auto.key=TRUE)
    
    # gráfico com legenda em colunas e alteração dos padrões dos pontos
    xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
           auto.key=list(columns=3),
           par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))
    
    # muda a posição da legenda
    xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
           auto.key=list(space="right", columns=1),
           par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))
    
    # coloca a legenda dentro do gráfico
    xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
           auto.key=list(corner=c(0.95,0.05), columns=1),
           par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))
    
    #-----------------------------------------------------------------------------
    # legenda para gráfico de linhas
    
    db <- with(da, aggregate(indice, list(cultivar=cultivar, dose=dose), mean))
    col <- trellis.par.get("superpose.line")$col[1:nlevels(db$cultivar)]
    
    # gráfico com legendas de linhas, nomes *antes* das linhas
    xyplot(x~dose, groups=cultivar, data=db, type="l",
           key=list(corner=c(0.95,0.05), columns=1,
             text=list(levels(db$cultivar)),
             lines=list(col=col)))
    
    # gráfico com legendas de linhas, nomes *depois* das linhas
    xyplot(x~dose, groups=cultivar, data=db, type="l",
           key=list(corner=c(0.95,0.05), columns=1,
             lines=list(col=col),
             text=list(levels(db$cultivar))))
    
    #-----------------------------------------------------------------------------
    # legenda para gráfico de linhas e pontos
    
    xyplot(x~dose, groups=cultivar, data=db, type="o", pch=19,
           key=list(corner=c(0.95,0.05), columns=1,
             type="o", divide=1,
             lines=list(col=col, pch=19),
             text=list(levels(db$cultivar))))
    
    #-----------------------------------------------------------------------------
    # legenda para pontos e reta de regressão
    
    # ajuste de modelo e criação de objeto com valores observados e preditos
    m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
    dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
    dc$indice <- predict(m0, newdata=dc)
    dc$ind <- "predito"
    da$ind <- "obsevado"
    dc <- rbind(da[,-3], dc)
    str(dc)
    
    # legenda fora do gráfico
    xyplot(indice~dose|cultivar, groups=ind, data=dc,
           distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
           key=list(space="top", columns=2, type="o", divide=1,
             lines=list(pch=c(1,NA), lty=c(0,1)),
             text=list(c("valores observados","valores preditos"))))
    
    # legenda dentro do gráfico
    xyplot(indice~dose|cultivar, groups=ind, data=dc,
           distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
           key=list(corner=c(0.95,0.05), columns=1, type="o", divide=1,
             lines=list(pch=c(1,NA), lty=c(0,1)),
             text=list(c("valores observados","valores preditos"))))
    
    #-----------------------------------------------------------------------------
    # dados de mortalidade de aves em função do tipo de resfriamento do aviário
    
    ma <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/mortes.txt",
                     header=TRUE, sep="\t")
    str(ma)
    
    #-----------------------------------------------------------------------------
    # legenda para gráficos de barras
    
    mb <- with(ma, aggregate(mortes, list(galpao=galpao, sistema=asper), sum))
    str(mb)
    
    col <- c("forestgreen", "palegreen")
    
    # legenda fora da região gráfica
    barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
             key=list(space="top", columns=2,
               rectangles=list(col=col),
               text=list(c("com aspersão","sem aspersão"))))
    
    # legenda fora da região gráfica com título
    barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
             key=list(title="Sistema de resfriamento", cex.title=1.1,
               space="top", columns=2,
               rectangles=list(col=col),
               text=list(c("com aspersão","sem aspersão"))))
    
    # legenda dentro da região gráfica com título
    barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
             key=list(title="Sistema de resfriamento", cex.title=1.1,
               corner=c(0.05,0.95), columns=1, padding.text=2,
               rectangles=list(col=col, size=3, height=0.8, border=TRUE),
               text=list(c("com aspersão","sem aspersão"))))
    
    #-----------------------------------------------------------------------------
    # escolhendo a posição da legenda com o mouse
    
    barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col)
    trellis.focus(name="page")
    loc <- grid.locator(unit="npc") # clicar no gráfico
    trellis.unfocus()
    
    barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
             key=list(title="Sistema de\nresfriamento", cex.title=1.1,
               x=loc$x, y=loc$y, corner=c(0.5,0.5), columns=1, padding.text=2,
               rectangles=list(col=col, size=3, height=0.8, border=TRUE),
               text=list(c("com\naspersão","sem\naspersão"))))
    
    #-----------------------------------------------------------------------------
    # nota de rodapé no gráfico com assinatura, data e versão
    
    # ajuste de modelo e criação de objeto com valores observados, preditos e IC
    m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
    dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
    indice <- predict(m0, newdata=dc, interval="confidence")
    dc <- cbind(dc, stack(data.frame(indice)))
    names(dc)[3] <- "indice"
    dc <- rbind(da[,-3], dc)
    
    #png("f017.png", w=500, h=350)
    xyplot(indice~dose|cultivar, groups=ind, data=dc,
           distribute.type=TRUE, type=c("l","l","p","l"),
           col=1, lty=c(1,2,0,2), layout=c(3,1),
           strip=strip.custom(bg="gray90"),
           xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
           key=list(space="top", columns=3, type="o", divide=1,
             lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
             text=list(c("valores observados","valores preditos","IC 95%"))))
    pushViewport(viewport())
    grid.text(label=paste("Walmes Zeviani --",
                format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
                R.version.string), rot=90,
              x=unit(1, "npc")-unit(2, "mm"),
              y=unit(0, "npc")+unit(2, "mm"),
              just=c("left", "bottom"),
              gp=gpar(cex=0.8, col="gray50"))
    popViewport()
    #dev.off()
    
    #-----------------------------------------------------------------------------
    # legenda para gráficos de 3 dimensões
    
    x <- seq(-10,10,l=50)
    dd <- expand.grid(x=x, y=x)
    dd$z <- with(dd, 500+x+2*y-0.5*x*y-x^2-y^2)
    
    # gráfico sem cores e sem legenda
    wireframe(z~x+y, data=dd)
    
    # gráfico com cores e legenda
    wireframe(z~x+y, data=dd, drape=TRUE)
    
    # remove a legenda
    wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE)
    
    # coloca a legenda embaixo
    wireframe(z~x+y, data=dd, drape=TRUE, colorkey=list(space="bottom"))
    
    # redimensionando a legenda
    wireframe(z~x+y, data=dd, drape=TRUE,
              colorkey=list(width=1, height=0.7))
    
    # gráfico de níveis de cor
    levelplot(z~x+y, data=dd)
    
    # com contornos e redimensiomento da legenda
    levelplot(z~x+y, data=dd, contour=TRUE,
              colorkey=list(width=1, height=0.7))
    
    #-----------------------------------------------------------------------------
    # usando a draw.key(), usar fora da função gráfica ou dentro da panel()
    
    xyplot(indice~dose|cultivar, groups=ind, data=dc,
           distribute.type=TRUE, type=c("l","l","p","l"),
           col=1, lty=c(1,2,0,2), layout=c(2,2),
           strip=strip.custom(bg="gray90"),
           xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico")
    draw.key(list(space="top", columns=1, type="o", divide=1,
             lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
             text=list(c("valores observados","valores preditos","IC 95%"))),
             draw=TRUE,
             vp=viewport(x=unit(0.75, "npc"), y=unit(0.65, "npc")))
    
    #-----------------------------------------------------------------------------
    # usando a draw.colorkey(), usar fora da função gráfica ou dentro da panel()
    
    wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE,
              panel=function(...){
                panel.wireframe(...)
                draw.colorkey(list(space="bottom", at=seq(250,500,25),
                                   height=0.8, width=1), draw=TRUE,
                              vp=viewport(x=unit(0.5,"npc"), y=unit(0.07,"npc")))
              })
    
    #-----------------------------------------------------------------------------
    # usando duas legendas
    
    de <- expand.grid(w=gl(3,5,la="w"), z=gl(4,3,la="z"))
    de$x <- with(de, rnorm(w, as.numeric(w), 0.5))
    de$y <- with(de, rnorm(z, as.numeric(z), 0.5))
    pch <- c(1,2,5,7); col=1:3
    
    xy <- xyplot(y~x, data=de, groups=z:w,
                 pch=rep(pch, e=3), col=rep(col, 4))
    
    print(xy, position=c(0,0,0.8,1), more=FALSE)
    draw.key(list(text=list(levels(de$z)), points=list(pch=pch, col=1),
             space="top", columns=1, title="Níveis de Z", cex.title=1.1),
             draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.6, "npc")))
    draw.key(list(text=list(levels(de$w)), points=list(pch=15, col=col),
             space="top", columns=1, title="Níveis de W", cex.title=1.1),
             draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.4, "npc")))
    
    #-----------------------------------------------------------------------------