Archive

Archive for the ‘delineamento’ Category

Fatorial duplo com dois tratamentos adicionais

Número de sementes em 40 infectadas por Fusarium (tranformado) em função da aplicação de fungicida em combinação com polímero.

Experimentos fatoriais são aqueles em que dois ou mais fatores estão presentes. Nos fatoriais incompletos nem todas as combinações possíveis estão presentes. Algumas das combinações podem não ser de interesse e por isso não são realizadas, alguma combinação é acidentalmente perdida durante o experimento ou certas combinações não identificam novos níveis, como por exemplo aplicação do produto A e B na dose 0 que, apesar dos rótulos diferentes, são o mesmo tratamento.

Um caso comum de fatorial incompleto são os fatoriais (completos) com a presença de tratamentos adicionais. Experimentos fatorais são realizados para verificar se os fatores interagem e então escolher combinação de níveis mais promissora para o processo estudado. Caso não exista interação, entende-se que os fatores têm efeito aditivos. Em qualquer caso, adicionam-se tratamentos adicionais ao experimento quando se deseja comparar o promissor aos de referência ou extremos.

Esse é um tipo de ensaio típico na avaliação do efeito de insumos agrícolas (fungicidas, insecticidas e herbicidas). Os ensaios testam produtos e suas doses e precisam ter as chamadas testemunhas para avaliar a eficiência das combinações geradas. Devido não ser um delineamento de níveis completamente cruzados a análise do experimento fatorial com tratamentos adicionais têm alguns detalhes que exigem mais atenção.

Nessa matéria faço análise de um fatorial com dois tratamentos adicionais. Os dados analisados foram retirados do Zimmermann (2004) de um experimento que estudou três fungicidas (benlate, captam, derosal) em três formas de aplicação (puro, antes de polímero, misturado com polímero). Adicionou-se dois tratamentos que foram a testemunha e aplicação do polímero puro. A resposta
avaliada foi o número de sementes infectadas com Fusarium em um total de 40 sementes.

O código está comentado para informar o leitor. São realizados a importação e reordenação dos níveis, a análise de variância com decomposição das somas de quadrados, contrastes entre efeitos, e testes sobre as médias ajustadas. Foi dada prioridade para os procedimentos matriciais de análise que podem ser adaptados para situações de desbalanceamento. Até a próxima ridícula.

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

require(lattice)
require(latticeExtra)
require(grid)
require(MASS)
require(doBy)
require(gmodels)
require(multcomp)

#-----------------------------------------------------------------------------
# leitura

fu <- read.table(
  "http://www.leg.ufpr.br/~walmes/data/zimmermann_fusarium.txt",
                 header=TRUE, sep="\t")
str(fu)

#-----------------------------------------------------------------------------
# modificar ordem dos níveis, primeiro níveis fatoriais, depois adicionais

levels(fu$fungicida)               # ordem correta
levels(fu$aplicacao)               # precisa alterar
levels(fu$aplicacao)[c(4,1,2,3,5)] # para ficar assim

fu$aplicacao <- factor(fu$aplicacao,
                       levels=levels(fu$aplicacao)[c(4,1,2,3,5)])
levels(fu$aplicacao) # agora ordem correta

#-----------------------------------------------------------------------------
# análise exploratória

xtabs(~fungicida+aplicacao, fu)

# os níveis não se cruzam completamente, portanto é um fatorial incompleto
# no caso 3x3+2, fatorial com 2 tratamentos adicionais

xyplot(infec40~tratamentos, fu) # resposta do tipo contagem de sucessos
xyplot(infec40~tratamentos, fu, jitter.x=TRUE)

#-----------------------------------------------------------------------------
# Vamos primeito analisar fora da estrutura fatorial, ou seja, considerando
# 3x3+2 = 11 níveis de um único fator.
# Vamos considerar distribuição para reproduzir os resultados do livro
# do Zimmermann (2004). Todavia, não é a análise mais correta devido os dados
# serem discretos. A distribuição binomial parece ser mais indicada.

m0 <- lm(infec40~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m0); layout(1) # não tão ruim, tentar boxcox?

m0 <- lm(infec40+1~tratamentos, data=fu) # soma 1 para eliminar os 0
boxcox(m0) # aponta log

m1 <- lm(log(infec40+1)~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m1); layout(1) # alguma melhoria?

# Zimmermann usa asin(sqrt(x/40)).
m2 <- lm(asin(sqrt(infec40/40))~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m2); layout(1) # alguma melhoria?

#-----------------------------------------------------------------------------
# Visualmente todos os resíduos são semelhantes, apresentam fugas.

lapply(list(m0, m1, m2), function(i) shapiro.test(residuals(i)))
lapply(list(m0, m1, m2), function(i) ks.test(scale(residuals(i)), "pnorm"))

# Pelos p-valores, a tranformação Box-Cox foi melhor, embora os resíduos ainda
# não apresentem um comportamento satísfatório.
# O foco não é testar transformações, então vamos seguir o que Zimmermann fez
# e trabalhar com asin(sqrt(x/40)).

#-----------------------------------------------------------------------------
# análise exploratória com a variável transformada

fu$y <- asin(sqrt(fu$infec40/40)) # cria a transformada

xyplot(y~tratamentos, fu, jitter.x=TRUE)

xyplot(y~fungicida, groups=aplicacao, fu, jitter.x=TRUE,
       type=c("p","a"), auto.key=TRUE)

# Apresenta indicativos de efeito de interação e que os adicionais
# são piores que o os fatoriais.

#-----------------------------------------------------------------------------
# análise na estrutura fatorial 3x3+2

m3 <- lm(y~fungicida*aplicacao, data=fu)
anova(m3)

# Os graus de liberdade são o que são pois:
# fungicida 4: são 5 niveis (3 fatoriais, 2 adicionais), 5-1=4;
# aplicacao 4: são 5, mas 2 estão confundidos (os adicionais), 5-2-1=2;
# fun:apli  4: (3-1)*(3-1)=4;

# O que queremos é a anova com somas de quadrados particionadas da seguinte
# forma
# fonte de variação     gl
# fungicidas             2
# aplicacao              2
# fun:apli               4
# fatorial vs adicionais 1
# polimero vs testemunha 1

#-----------------------------------------------------------------------------
# Na minha opinião, obter esse quadro de análise de variância é desperdício
# de tempo pois:
# * não tem sentido prático comparar o efeito médio dos fatoriais contra o
#   efeito médio dos adicionais, principalmente sabendo que a interação
#   fungicida:aplicacao foi significativa;
# * a interação fungicida:aplicacao significativa torna não interpretável
#   o teste para os efeitos principais;
# * polímero vs testemunha tem só um grau de liberdade, posso comprar os
#   efeitos por um testes t;

coef(m3) # tem NA para os efeitos não estimados (não presentes)
efstm <- colnames(vcov(m3)) # nomes dos efeitos que foram estimados
b <- coef(m3)[efstm] # só os estimados
b

# matriz de coeficientes para estimação das médias ajustadas
X <- popMatrix(m3, effect=c("fungicida", "aplicacao"))
str(X)

Xn <- attr(X, "grid") # 25 linhas, 5x5=25 (faz pensando que é completo)
Xo <- unique(fu[,c("fungicida","aplicacao")]) # que ocorrem (3x3+2)
Xo

#-----------------------------------------------------------------------------
# calculando as médias ajustadas matricialmente

rownames(X) <- apply(Xn, 1, paste, collapse=":")
on <- apply(Xo, 1, paste, collapse=":")
X <- X[on,efstm] # só linhas e colunas que correspondem à efeitos que ocorrem

X%*%b                                             # são as médias ajustadas
with(fu,
     tapply(y, list(fungicida, aplicacao), mean)) # são médias amostrais

# elas coincidem por ser um experimento regular (balanceado/ortogonal)

#-----------------------------------------------------------------------------
# fazendo contrastes
# não funcionam com presenta de NA
# contrast::contrast()
# multcomp::glht()
# para isso usar a gmodels::estimable()

m4 <- aov(formula(m3), data=fu) # aov não NA em coef

# polímero contra a testemunha
p.vs.t <- matrix(X[10,]-X[11,], nrow=1)
rownames(p.vs.t) <- "p.vs.t"
estimable(m4, cm=p.vs.t)

# fatorial contra adicional
f.vs.a <- matrix(apply(X[1:9,], 2, mean)-apply(X[10:11,], 2, mean), 1)
rownames(f.vs.a) <- "f.vs.a"
estimable(m4, cm=f.vs.a)

estimable(m4, cm=rbind(p.vs.t, f.vs.a)) # os dois de uma só vez

# E aí, complicado fazer tudo isso, né? Qual a sensação em saber de que
# pouco disso vai ser útil? Se deu interação teremos que estudar níveis
# de um fator fixando os níveis do outro.

#-----------------------------------------------------------------------------
# fazendo a decomposição das somas de quadrados dentro da anova

cf <- contrasts(C(fu$fungicida, helmert))
cf # o segredo está em montar os contrastes corretamente
cf[,3] <- c(-1,-1,-1,1.5,1.5) # fatorial contra adicional
cf[,4] <- c(0,0,0,-1,1)       # polímero vs testemunha
cf
contrasts(fu$fungicida) <- cf

m5 <- aov(y~fungicida*aplicacao, data=fu)
summary(m5, expand=FALSE,
        split=list("fungicida"=list(
                     "fungicida"=1:2,
                     "fat.vs.adi"=3,
                     "pol.vs.tes"=4)))

# Esse procedimento depende de balanceamento pois as somas de quadrados são
# sequênciais.

#-----------------------------------------------------------------------------
# Uma vez que deu interação é de prática fazer o desdobramento que consiste
# em estudar o efeito de aplicação dentro dos níveis de fungicida e
# vice-versa. Além do mais, por ter tratamentos adicionais, pode-se comparar
# com os níveis dos fatoriais. Vamos listas o número de hipóteses que teremos
# que fazer:
# * 3 níveis de aplicação em cada fungicida: 3*choose(3,2)=9
# * 3 níveis de fungicida em cada aplicação: 3*choose(3,2)=9
# * testemunha contra cada um dos 9 níveis do 3x3: 9
# * polimero contra cada um dos 9 níveis do 3x3: 9
# Com isso são 4*9=36 comparações. Já que é assim, vou fazer todos
# contra todos, choose(11, 2)=55, que o usuário pode olhar a que lhe
# interessa, afinal de contas, de 36 para 55 não têm diferença

mm <- model.matrix(m4)
dim(mm)
mm <- mm[,efstm]

# modelo de classe lm equivalente à m5, sem NA
m6 <- lm(y~-1+mm, data=fu)
anova(m6)

cbind(coef(m4), coef(m6)) # são o mesmo modelo

#-----------------------------------------------------------------------------
# médias ajustadas (ma)

cima <- confint(glht(m6, linfct=X)) # intervalos de confiança para ma
cima # tem cobertura *global* de 95% (e não indivídual)

str(cima)
cima$confint

result <- cbind(Xo, cima$confint)
str(result)

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm")

#-----------------------------------------------------------------------------
# carrega algumas funções que serão úteis para comparar médias e fazer
# gráficos com intervalos de confiança

objs <- ls()
objs <- c(objs, "apc", "prepanel.cbH", "panel.cbH")
source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functionsi.R",
       encoding="latin1")

rm(list=ls()[!(ls()%in%objs)])
ls()

#-----------------------------------------------------------------------------
# comparando as médias ajustadas, choose(11, 2), contrastes de Tukey

str(X)
Xc <- apc(X)
str(Xc)

#-----------------------------------------------------------------------------
# fazendo as comparações com a glht, método para correção de p-valor
# fdr, o single-step demora muito

c0 <- summary(glht(m6, linfct=Xc), test=adjusted(type="fdr"))
c0

c0$focus <- "comparacoes" # para poder usar a cld()
cld(c0)

result$cld <- cld(c0)$mcletters$Letters

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

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm",
        panel=function(x, y, z, centers, subscripts, ...){
          panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...)
          panel.text(centers, z,
                     paste(format(centers, digits=2),
                           result$cld[subscripts], sep=" "), pos=3)
        })

#-----------------------------------------------------------------------------
# gráfico final

result$desloc <- 0.05*c(0,-1,1,0,0)[match(result$aplicacao,
                                          levels(fu$aplicacao))]

#png("f048.png", 500, 400)
xyplot(y~fungicida, groups=aplicacao, data=fu, jitter.x=TRUE, amount=0.05,
       xlab="Fungicida", ylab=expression(asin(sqrt(x/40))),
       prepanel=prepanel.cbH, ly=result$lwr, uy=result$upr,
       auto.key=list(title="Aplicação", cex.title=1.2, columns=2,
         type="o", divide=1, lines=TRUE, points=FALSE))+
  as.layer(xyplot(Estimate~fungicida, groups=aplicacao, data=result, type="l",
                  ly=result$lwr, uy=result$upr,
                  desloc=result$desloc,
                  cty="bars",
                  panel.groups=panel.cbH,
                  panel=panel.superpose))
trellis.focus("panel", 1, 1, h=FALSE)
x <- as.numeric(result$fungicida)+result$desloc
y <- result$Estimate
lab <- paste(format(y, digits=2), result$cld, sep=" ")
a <- paste(rep("0", max(nchar(lab))), collapse="")
grid.rect(x=unit(x, "native"), y=unit(y, "native"),
          width=unit(1, data=a, "strwidth"),
          height=unit(1.2, data=a, "strheight"),
          hjust=-0.1, gp=gpar(fill="gray90", col=NA, fontsize=10))
grid.text(lab, x=unit(x, "native"), y=unit(y, "native"),
          hjust=-0.2, gp=gpar(col="black", fontsize=10))
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
Categorias:delineamento Tags:, ,

Semântica para descrever modelos

A função lm() do R resolve 90% dos casos de análise de experimentos. Não são para lm() modelos não lineares, de efeito aleatório, de distribuição diferente da gaussiana ou dados não independentes (correlacionados). Para uso eficiente da função é importante saber declarar os modelos. Para isso existe uma semântica (termo bem apropriado usado pelo Fernando Toledo) para declarar os modelos. Aqui eu vou listar alguns dos casos possíveis.

fórmula significado
A+B efeito simples aditivo entre fatores
A:B efeito do produto cartesiano de níveis dos fatores
A*B*C efeitos simples e todas interações possíveis
(A+B+C)^2 efeitos simples e interações de até 2 termos (duplas)
(A+B+C+C+E)^3 efeitos simples e interações de até 3 termos (triplas)
A/B efeito simples de A e de B aninhado em A, A+A:B
I(x^2) operador I() torna x^2 literal (remove o poder semântico de ^2)
-A:B remove o termo do modelo
A*B*C-A:B:C declara modelo triplo mas remove interação tripla, (A+B+C)^2
0+A remove intercepto (coluna de 1) do modelo
-1+A idem ao anterior
poly(x, degree=2) polinômio ortogonal de grau 2 em x
poly(x, degree=2, raw=TRUE) polinômio cru de grau 2 em x
A*(B+C) propriedade distributiva, A*B+A*C
A/(B+C) idem ao anterior, A/B+A/C
A/B/C propriedade recursiva, A+A:B+A:B:C
. representa todas as colunas, exceto a da resposta
(.)^2 inclui todas as interações duplas das variáveis

Além dessas opções de sintaxe, ainda temos a possibilidade de usar o operador semântico Error() para declarar termos de efeito aleatório. Esse operador só é interpretado pela função aov(). Com ele fazemos a análise de experimentos em parcelas subdivididas, por exemplo, de tal forma que a análise de variância apresenta os testes F considerando os quadrados médios corretos.
Além disso podermos usar a função terms(..., keep.order=TRUE) para fixar a ordem dos termos na análise de variância, pois os termos são ordenados pela ordem de magnitude, assim, sempre interações vem depois de efeitos simples, algumas vezes não desejamos isso, como é o caso da análise de experimentos em blocos incompletos. Os exemplos vão tornar esses apontamentos mais claros. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# observe as matrizes de delineamento geradas nessas situações

da <- replicate(5, gl(3,1), simplify=FALSE)
da <- data.frame(do.call(expand.grid, da))
names(da) <- LETTERS[1:ncol(da)]
da$x <- runif(nrow(da))

formulas <- list(~A+B, ~A:B, ~A*B, ~(A+B)^2, ~(A+B+C+D+E)^3,
                 ~A/B, ~A/B/C, ~A/(B+C), ~A*(B+C), ~A*(B+C)^2,
                 ~(A+B+C)^2-A:B, ~x, ~0+A, ~-1+A, ~0+x, ~A*x, ~A:x,
                 ~A/x, ~x+I(x^2), ~poly(x, degree=2),
                 ~poly(x, degree=2, raw=TRUE), ~A+A:B+C,
                 terms(~A+A:B+C, keep.order=TRUE))
names(formulas) <- apply(sapply(formulas, paste), 2, paste, collapse="")

X <- lapply(formulas, model.matrix, data=da)
lapply(X, head)

#-----------------------------------------------------------------------------
# para obter a equação de um polinômio para os níveis

soja <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)

soja$A <- factor(soja$agua)
soja$K <- soja$potassio

m0 <- lm(rengrao~bloco+A*(K+I(K^2)), soja) # I() torna literal
m1 <- lm(rengrao~bloco+A*poly(K, degree=2, raw=TRUE), soja)  # original
m2 <- lm(rengrao~bloco+A*poly(K, degree=2, raw=FALSE), soja) # ortogonal

cbind(coef(m0), coef(m1), coef(m2)) # 2 primeiras colunas são iguais

m3 <- lm(rengrao~0+A/(K+I(K^2))+bloco, soja) # para ter os coeficientes
coef(m3) # os interceptos, termos lineares e termos quadráticos
c0 <- coef(m3)[grep("bloco", names(coef(m3)), invert=TRUE)]
matrix(c0, ncol=3, dimnames=list(levels(soja$A), c("b0","b1","b2")))

#-----------------------------------------------------------------------------
# tranformações dentro da fórmula

m0 <- lm(rengrao~bloco+A*(K+sqrt(K)), soja)  # não requer I()
m0 <- lm(rengrao~bloco+A*I(K-mean(K)), soja) # centrar na média

#-----------------------------------------------------------------------------
# fazendo de conta que o experimento é em parcelas subdividas para ilustrar
# o uso do operador semântico Error()

soja$K <- factor(soja$K)

# combinar os níveis de bloco e A nos identificamos as parcelas
m4 <- aov(rengrao~bloco+A*K+Error(bloco:A), soja) # dá mensagem de aviso
summary(m4) # anova separa pelos estratos

soja <- transform(soja, parcela=interaction(bloco, A))

# é o mesmo modelo
m5 <- aov(rengrao~bloco+A*K+Error(parcela), soja) # não dá mensagem de aviso
summary(m5) # anova separa pelos estratos

#-----------------------------------------------------------------------------
# fazendo uso da função terms() para fixar a ordem dos termos na anova()

bib1 <- read.table("http://www.leg.ufpr.br/~walmes/data/bib1.txt",
                   header=TRUE, sep="\t")
str(bib1)

bib1 <- transform(bib1,
                  repetição=factor(repetição),
                  variedade=factor(variedade),
                  bloco=factor(bloco))
str(bib1)

# repetição + bloco dentro de repetição + variedade
m0 <- lm(y~repetição/bloco+variedade, bib1)
anova(m0) # por ser de maior ordem o termo repetição:bloco é o último

m1 <- lm(terms(y~repetição/bloco+variedade, keep.order=TRUE), bib1)
anova(m1) # ordem obedeceu a ordem dentro da fórmula

#-----------------------------------------------------------------------------
Categorias:delineamento Tags:, ,

PCA de três dimensões

Olá a todos. O objetivo desse post é mostrar como fazer a análise de componentes principais em três dimensões. Mas esse post tem o destaque não ter sido escrito sozinho… Sem a participação da Msc Thalita (ou “Hermínia”, ninguém vai entender!) não teria post. Agradeçam a ela, mas reclamem comigo! Pois bem, Vamos usar como exemplo dados de um estudo com fatorial duplo (4 \times 4), em que foram tomadas quatro características (dados disponíveis para teste!).

Vou aproveitar e deixar aqui também um alô para o Dsc Guilherme pela ajuda nos detalhes do script, Valeu!

Um pouco de teoria:

Experimentos fatoriais (n \times n) caracterizam-se pela apresentação dos dados em tabelas de duas entradas (matriz), sendo que cada casela da tabela contém a resposta média de cada combinação dos fatores sob análise. Todas possibilidades desses experimentos estão fora do escopo desse post, por isso vamos um pouco além…

Acrescentando-se um outro fator, ou ainda outra característica (n \times n \times n, por exemplo), com 3 interações duplas e uma tripla, isso fica um pouco mais complexo. Neste caso, os dados são organizados em “CUBO” (arranjo de três entradas).

Nesta situação não é possível aplicar a análise de componentes principais usual. Uma alternativa é a utilização da decomposição PCAn, proposta por Tuckey em 1964.

Portanto nesse post pretendemos (afinal esse foi um trabalho em equipe) mostrar como realizar a decomposição de um arranjo de três entradas (CUBO). E mais que isso, mostrar uma função que realiza a escolha entre os vários modelos possíveis que variam conforme a dimensão dos dados.

Vamos ao código… mas calma, primeiro preparamos os dados!

Veja que usamos um pacote (PTAk, disponível no CRAN). Caso você não o tenha, proceda com install.packages(‘PTAk’)!

require(PTAk) # pacote analises pscicrometricas

## le e transforma os dados em fatores
dados <- transform(read.table("http://dl.dropbox.com/u/38195533/dados.txt")
                              na.string = '.',
                              header = TRUE,
                              dec = ',',
                              sep = '\t'),
                   elemento = factor(elemento),
                   residuo = factor(residuo),
                   metodo = factor(metodo))

CUBO.inicial <- with(dados,
                     tapply(conc, list(residuo, elemento, metodo),
                            mean, na.rm = TRUE))
## padronização dos dados
CUBO.standard <- apply(CUBO.inicial, 2, scale, center = FALSE)

## reconstrucao da planilha
estrutura <- expand.grid(elemento = factor(1:4),
                         residuo = factor(1:4),
                         metodo = factor(1:4))

plan <- data.frame(estrutura, conc = c(CUBO.standard)) ## dados

ajuste <- aov(conc ~ elemento * residuo * metodo - elemento:residuo:metodo,
              data = plan)
erros <- ajuste$residuals

## CUBO
Z <- tapply(erros, list(plan$residuo, plan$elemento, plan$metodo), mean) 

Nesse trecho são lidos os dados, na forma usual de uma planilha de dados (uma coluna para cada variável). Pela natureza dos dados, estes são padronizados (N(\mu = 0, \sigma^2 = 1)). Depois a panilha de dados é “reconstruída” e feito o ajuste do modelo linear respectivo (sem a interação tripla — objeto do nosso estudo!). Do resultado da análise extraímos os resíduos e construímos o CUBO.

Quase tudo pronto! O próximo trecho é a função que efetivamente faz os vários ajustes possíveis (isso é função da dimensão dos dados).

otimim.PCAn <- function(cubo) {
  dimensoes <- dim(Z) - 1
  diagnostico <- vector('list', length = prod(dimensoes - 1))
  ajustes <- vector('list', length = prod(dimensoes - 1))
  contador <- 0
  for(i in 2:dimensoes[ 1 ])
    {
      for(j in 2:dimensoes[ 2 ])
        {
          for(k in 2:dimensoes[ 3 ])
            {
              modelo <- c(i, j, k)
              tamanho <- sum(modelo)
              ## decomposicao Tucker -- PTAk
              tucker <- PCAn(Z, # array
                             dim = modelo, # dimensoes
                             test = 1E-12,
                             Maxiter = 1000,
                             smoothing = FALSE,
                             verbose = FALSE,
                             file = NULL)
              ## diagnostico
              porcentagem <- tucker[[ 3 ]]$pct # % explicada
              contador <- contador + 1
              ajustes[[ contador ]] <- tucker
              diagnostico[[ contador ]] <- c(modelo, tamanho, porcentagem)
            }
        }
    }
  ## diagnosticos dos ajustes
  tab.ajustes <- data.frame(do.call(rbind, diagnostico))
  tab.ajustes <- tab.ajustes[order(-tab.ajustes[, 4], -tab.ajustes[, 5]), ]
  por.total <- split(tab.ajustes, tab.ajustes[, 4])
  melhores <- vector('list', length = length(por.total))
  for(i in 1:length(melhores))
    {
      melhores[[i]] <- por.total[[i]][1, ]
    }
  ## melhores ajustes por dimensao
  tab.melhores <- do.call(rbind, melhores)
  tab.melhores <- tab.melhores[order(-tab.melhores[, 4]), ]
  ## melhor modelo
  melhor <- tab.melhores[(which((tab.melhores[,5]-
                          c(tab.melhores[,5][-1],0))>=5)[1]+1),]
  ## resposta -- melhor ajuste
  return(list(escolhido = melhor,
              resposta = ajustes[[ as.numeric(rownames(melhor)) ]]))
}

UFA!

Feita a função agora é só usa-lá (simples não?). Um uso possível dessa decomposição é representar e interpretar os dados graficamente em um biplot. Nessa condição de arranjos de três entradas são possíveis diferentes formas de gráficos biplots. Mas o objetivo do nosso post não é apresentá-los.

O trecho final aplica a função recém-construída (e retorna um controle de tempo de cada ajuste). Mostramos também como extrair as matrizes (A, B e C) e o arranjo G do melhor ajuste.

As matrizes de resposta e o arranjo G é que são úteis na construção dos gráficos. Aproveitem!

ajuste.PCA <- otimim.PCAn(Z)

A <- ajuste.PCA$resposta[[1]]$v
B <- ajuste.PCA$resposta[[2]]$v
C <- ajuste.PCA$resposta[[3]]$v

G <- ajuste.PCA$resposta[[3]]$coremat

summary(ajuste) # ANOVA do modelo
summary.PTAk(ajuste.PCA$resposta) # quanto cada componente explica

Voi lá…, Até mais pessoal!

Experimento fatorial duplo com um tratamento adicional

Produção de cultivares em função do nível de fertilizante (95% de confiança). As linhas horizontais em cinza representam o desempenho da cultivar de referência. Linhas pretas representam o desempenho das demais como função do nível de fertilizante. A seta aponta o nível de fertilizante que confere a máxima produção da cultura. Análise baseada em dados simulados.

Em algumas áreas do conhecimento é quase regra fazer experimentos incluindo uma tratamento adicional. Um experimento fatorial é aquele que possui duas ou mais fontes de variação. Em geral, têm-se presente no experimento todos os níveis possíveis das combinações entre os níveis desses fatores, ou em outras palavras, todos os níveis possíveis do produto cartesiano. Os fatores podem ser qualitativos (ou categóricos) ou quantitativos (ou métricos). Um tratamento adicional é um nível que não corresponde à estrutura fatorial, ou seja, é um nível que não veio da combinação dos níveis dos fatores. Para ficar mais claro, segue alguns exemplos:

  • Para um determinado município, sabe-se que a melhor cultivar de soja é a A com produção máxima ao receber a aplicação de 100 kg de fertilizante. Uma empresa está estudando 5 novas cultivares para o município e precisa definir a dose de fertilizante para a máxima produção. Além de definir a dose ótima, deseja-se comparar a produção com a cultivar A. Para isso faz-se um experimento com os 5 níveis de cultivar de soja (B, C, D, E, F) e um número de níveis suficientes de fertilizante que permita determinar o ponto de máximo, diga-se, 5 níveis (50, 75, 100, 125, 150 kg). O experimento terá 25 níveis (5×5) vindos da porção fatorial e um nível de referência, chamado de testemunha ou controle, que é cultivar A com 100 kg de fertilizante. Esse experimento é um fatorial 5×5+1.
  • Em um experimento com herbicidas, deseja-se definir a dose mínima de aplicação para controlar plantas daninhas. São 4 níveis de novos herbicidas (A, B, C, D) e 6 níveis de concentração (0,1; 0,5; 1; 2; 5; 10%). A variável resposta de interesse é produtividade da cultura. Para verificar se os produtos causam toxidade à cultura, tem-se um nível que é capina manual (testemunha positiva), que é o controle de plantas daninhas sem produto químico, onde se espera a maior produtividade. Para ver a eficiência desses herbicidas, no estudo se incluí um nível que é sem a capina (testemunha negativa), onde se espera a menor produtividade. Este experimento é um fatorial 4×6+2.

Perceba que pelas descrições acima ficaram evidentes às hipóteses de cada experimento. No primeiro é definir a dose ótima de cada cultivar e depois comparar o desempenho das cultivares na melhor dose com a cultivar A. No segundo é verificar se os herbicidas causam toxides à cultura (comparar com a testemunha positiva), qual a perda de produção quando não se faz o controle de plantas daninhas (comparar com a testemunha negativa) e definir a dose mínima para controlar plantas daninhas para cada herbicida.

O modelo estatístico que representa o primeiro experimento é

Y_{ijk} = M+C_i + F_j + CF_{ij} + \epsilon_{ijk} \quad \text{e} \quad Y_{k} = M + A + \epsilon_{k},

em que Y são os valores observados, M um valor inerente à todas as observações, C representa o efeito dos níveis de cultivar, F representa o efeito dos níveis de fertilizante, CF representa o efeito da combinação dos níveis de cultivar e fertilizante, A é o efeito do tratamento adicional e \epsilon é um erro aleatório.

Veja que até aqui foi tudo muito comum, com exceção ao fato de ser necessário escrever o modelo com duas expressões, uma para descrever a parte fatorial e outra para a parte adicional. No R não conseguimos declarar o modelo com duas fórmulas e é justamente sobre isso que vou tratar nesse post.

Podemos representar o modelo acima como um fatorial incompleto. Pense assim: são 6 nível de cultivar e apenas 5 deles combinados com 5 níveis de fertilizante. É assim que vamos declarar o modelo no R. Com algumas manipulações, iremos repartir a soma de quadrados da porção fatorial da porção adicional. A repartição será ortogonal que nem sempre representa as hipóteses que nós estabelecemos. Então por que fazemos essa repartição? Só porque ela é ortogonal.

Vou usar dados simulados para exemplificar o procedimento. O efeito dos níveis de fertilizante serão representados por um polinômio de segundo grau, apenas para que apresente ponto de máximo. No script abaixo eu gero os dados, faço os ajuste dos modelos, apresento os quadros de análise de variância, faço a repartição da soma de quadrados, a estimação da dose ótima de fertilizante para a resposta máxima e a comparação com a testemunha. O código foi comentado para orientar o leitor. Se você quer saber como é a análise do fatorial com dois tratamentos adicionais, deixe um comentário nesse post. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# gerando dados conforme o modelo estatístico 5x5+1, depois explico o fert=101

adi <- expand.grid(cult=gl(1,5,la=LETTERS[1]), fert=101)
fat <- expand.grid(cult=gl(5,5,la=LETTERS[2:6]), fert=seq(50,150,25))
da <- rbind(adi, fat)

theta <- c(c(-194.29, -197.26, -197.85, -203.03, -190.20, -190.45),
           c(9.1797, 8.2686, 8.6437, 9.3438, 8.8773, 8.1872),
           c(-0.03382, -0.03479, -0.03632, -0.03341, -0.03597, -0.03675))
X <- model.matrix(~-1+cult/(fert+I(fert^2)), data=da)
da$eta <- X%*%matrix(theta)
set.seed(2)
da$y <- da$eta+rnorm(nrow(da),0,30)

require(lattice)
xyplot(y~fert|cult, data=da, type=c("p","a"))

#-----------------------------------------------------------------------------
# usar o confundimento

da$Fert <- factor(da$fert) # esse passo explica o fert=101
levels(da$Fert) # o nível 101 ocorre nas mesmas observações onde
levels(da$cult) # ocorre o nível A, portanto estão confundidos

#-----------------------------------------------------------------------------
# ajuste do modelo considerando os fatores todos de níveis categórios
# poderia ir direto para a análise com polinômio mas esse passo ilustra
# o procedimento para 2 fatores de níveis categóricos

m0 <- aov(y~cult*Fert, data=da)
anova(m0) # compare os graus de liberdade com o número de níveis

# 6-1 gl para cult
# 5-1 gl para Fert
# (6-1)*(5-1) para cult:Fert? não, é (5-1)*(5-1)
# pois existe níveis de fert em apenas 5 níveis de cultivar

#-----------------------------------------------------------------------------
# os contrates entre níveis são os seguintes

contrasts(da$cult) # todo nível de cult contra o nível A
contrasts(da$Fert) # todo nível de Fert contra o nível 50

#-----------------------------------------------------------------------------
# para decompor a soma de quadrados em parte fatorial e parte adicional
# devemos mudar esses contrastes. Um deles deve ser o nível adicional vs todos
# os demais níveis daquele fator. Para isso podemos trocar a ordem dos níveis
# e o tipo de contraste para Helmert

contrasts(C(da$cult, contr=helmert)) # nível F contra os demais, trocar para A
contrasts(C(da$Fert, contr=helmert)) # nível 150 trocar para 101

da$cult <- factor(da$cult, levels=rev(levels(da$cult)))
da$Fert <- factor(da$Fert, levels=c(50,75,100,125,150,101))

levels(da$cult) # A é o último nível
levels(da$Fert) # 101 é o último nível

contrasts(C(da$cult, contr=helmert)) # A vs demais
contrasts(C(da$Fert, contr=helmert)) # 101 vs demias

m0 <- aov(y~cult*Fert, data=da, # passando o tipo de contraste
          contrasts=list(cult=contr.helmert, Fert=contr.helmert))
anova(m0) # mesmo quadro de anova

# agora vamos repartir a soma de quadrados. Note que existe uma hipótese em
# teste quando fazemos esse procedimento. Essa hipótese surge porque estamos
# fazendo um partição ortogonal específica.

summary(m0, expand=FALSE, split=list("cult"=list(fat=1:4, adi=5)))

# de 1:4 são estimativas dos contrastes entre níveis da porção fatorial
# 5 é a estimativa do contraste entre o nível adicional contra os outros
# e essa é a hipótese que a partição ortogonal testa
# mas que nem sempre é a hipótese que temos interesse
# agora é só aplicar os testes para as hipótese de interesse, sejam elas
# comparações múltiplas, contrastes planejados, etc.

#-----------------------------------------------------------------------------
# ajustando modelo polinomio de 2 grau para o efeito de fertilizante

m1 <- lm(y~-1+cult/(fert+I(fert^2)), data=da)
summary(m1)

# usamos essa fórmula para termos interpretação direta das estimativas
# perceba que para o nível A obtivemos NA, lógico, não temos como estimar
# efeito de fertilizante porque não temos níveis para essa cultivar

#-----------------------------------------------------------------------------
# estimando as doses que conferem máxima resposta em cada cultivar
# solução é derivar o modelo quadrático e igualar a zero

D(expression(b0+b1*x+b2*x^2), "x") # igualar a zero e isolar x
fert.max <- -0.5*coef(m1)[6+1:6]/coef(m1)[12+1:6]
fert.max <- fert.max[-6]
fert.max # as doses de fertilizante que conferem o máximo para cada cultivar

#-----------------------------------------------------------------------------
# vamos ajustar com a aov() que omite os NA para então usar a estimable()

m2 <- aov(y~-1+cult/(fert+I(fert^2)), data=da)
summary.lm(m2) # não contém NA

#-----------------------------------------------------------------------------
# agora é só estimar a produção máxima para cada nível de cultivar passando a
# para estimable() a matrix que identifica as funções lineares envolvidas

con <- cbind(diag(5), 0, diag(5)*fert.max, diag(5)*fert.max^2)
con <- rbind(con, 0); con[6,6] <- 1
rownames(con) <- paste(levels(da$cult),
                       c(round(fert.max,1),"100.0"), sep="-")
con

require(gmodels)
e0 <- estimable(m2, cm=con, conf.int=0.95)
e0 # quadro de estimativas da produção nas doses de máximo

#-----------------------------------------------------------------------------
# fazer um gráfico da produção máxima com os intervalos

xlim <- extendrange(e0[,6:7])
barchart(rownames(e0)~Estimate, data=e0, xlim=xlim,
         xlab="Produção", ylab="Cultivar-dose",
         panel=function(x,y,...){
           panel.barchart(x,y,...)
           panel.arrows(e0$Upper.CI, y, e0$Lower.CI, y,
                        code=3, angle=90, length=0.1)
         })

#-----------------------------------------------------------------------------
# fazendo as comparações das produções máximas com a testemunha

contr <- t(sapply(1:5, function(l) con[l,]-con[6,]))
rownames(contr) <- paste(rownames(con)[-6], "vs A")
e1 <- estimable(m2, cm=contr, conf.int=0.95)
e1

# segundo os contrastes, apenas o cultivar D é superior ao A quando recebe
# a dose 158 de fertilizante

#-----------------------------------------------------------------------------
# confeccionando o gráfico final

pred <- expand.grid(cult=levels(da$cult)[-6], fert=seq(50,170,1))
aux <- predict(m2, newdata=pred, interval="confidence")
aux <- stack(as.data.frame(aux))
pred <- cbind(pred, aux)
str(pred)

require(latticeExtra)

#png("f020.png", w=600, h=300)
xyplot(y~fert|cult, subset(da, cult!="A"), col="green4",
       xlim=extendrange(pred$fert), ylim=c(0,525), layout=c(5,1),
       xlab="Nível de fertilizante", ylab="Produtividade da cultura",
       strip=strip.custom(bg="gray90"),
       panel=function(...){
         panel.xyplot(...)
         panel.abline(h=e0[6,c(1,6,7)], lty=c(1,2,2), col="gray50")
         panel.arrows(fert.max[which.packet()], e0[which.packet(),1],
                      fert.max[which.packet()], 0, length=0.1)
       })+
  as.layer(xyplot(values~fert|cult, groups=ind, data=pred,
                  type="l", col=1, lty=c(1,2,2)))
trellis.focus("panel", 1, 1)
panel.text(50, e0[6,1], label="Cultivar A", pos=4, font=3, cex=0.8)
trellis.unfocus()
#dev.off()

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

Perda de observações em experimento em DBC

Muitas pessoas não sabem o que fazer quando tem observações não disponíveis no seu experimento – o termo observações perdidas dá a sensação de que um dia elas existiram mas por descuido se perderam. Outras pessoas sentem-se seguras porque têm a licença de determinado aplicativo para fazer análise que, segundo a crença popular, resolve tudo. Bem, o que se deve ter claro em mente, em qualquer condição experimental é: quais as hipóteses de interesse e quais os pressupostos.
Um experimento em delineamento de blocos ao acaso (os blocos não são sorteados, mas os níveis dos demais fatores que são sorteados em cada bloco) com um único fator, normalmente chamado fator tratamento, é representado pelo seguinte modelo estatístico

Y_{ij} = \mu + \beta_i + \tau_j + \epsilon_{ij} ,

em que \mu é uma constante incidente em todas as observações, \beta representa o efeito dos níveis do fator bloco, \tau representa o efeito dos níveis do fator tratamento e \epsilon é um termo que representa o erro. Os pressupostos desse modelo são:

  • os efeitos de bloco e tratamento são aditivos;
  • o erro é independente, de mesma variância e com distribuição normal (condição necessária para condução de inferências);

O fator bloco representa alguma fonte de variação identificável nas unidades experimentais (UE), como por exemplo idade, tamanho, vigor, topografia, umidade, procedência, posição, histórico, comportamento, sexo. Caso não exista uma fonte de variação identificável não há porque, nem como classificar/separar, as UE em blocos que sejam diferentes entre si. Há casos em que as UE são iguais mas ao longo do período experimental podem sofrer influências identificáveis devido, por exemplo, às pessoas que operam, as posições das UI no campo, às épocas de colheita, aos avaliadores. As diferenças entre as UE de um mesmo bloco não podem ser causadas por uma fonte identificável, e sim serem devidas somente ao acaso. Em cada nível de bloco, os níveis de tratamento são sorteados às UE. O experimento é conduzido e as variáveis resposta são observadas.

Às vezes, ocorre de algumas variáveis não serem medidas em algumas UE. Se os valores não disponíveis (NA, not available) surgiram ao acaso, a análise não é inviabilizada. Há casos em que NA depende do tratamento, por exemplo, vacas invadem a área experimental e comem as UI que são do melhor capim, uma doença não controlada mata as UE do genótipo não resistente, pessoas desavisadas colhem os frutos da variedade precoce do pomar experimental. Nesse caso, os NA não sugiram ao acaso e a análise não será confiável.

Nesse tipo de experimento, não interesse em avaliar as diferenças entre níveis de bloco. Já espera-se haver diferenças mas não se quer fazer inferências sobre esses valores. Os blocos são (e devem ser) considerados na análise porque são uma fonte identificável de variação. Por outro lado, para os níveis do fator tratamento se estabelece um conjunto de hipóteses de interesse. São comuns comparações entre os valores esperados (\mu+\tau ) dos níveis contra o nível de referência (testemunha/controle) e comparações duas a duas. Então é necessário que se estime os parâmetros para testar essas hipóteses. No caso balanceado, a média amostral para os níveis de tratamento é um ótimo estimador para (\mu+\tau ) mas caso desbalanceado, a média amostral é um estimador viesado.

O procedimento que vou apresentar aqui é matricial. É o procedimento geral para estimação das médias populacionais marginais (ou valor esperado, \mu+\tau ). As médias amostrais coincidem com esse estimador matricial no caso balanceado. Esse estimador é chamado de estimador de mínimos quadrados. No final apresento um ajuste considerando os efeitos de bloco aleatório mas que vou detalhar as propriedades num post futuro. O código abaixo está comentado para orientação do leitor. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados artificiais gerados a partir do modelo de 4 blocos e 3 tratatamentos

da <- expand.grid(bl=gl(4,1,la="b"), tr=gl(3,1,la="t"))

# vetor de parâmetros (ou efeitos)
theta <- c(mu=10, b1=-1, b2=-1, b3=0, b4=2, t1=-1, t2=0, t3=1)

# soma dos efeitos de bl e tr dá zero, isso é apenas para facilitar
sum(theta[2:5]); sum(theta[6:8])

# médias populacionais marginais dos níveis do fator tratamento
# obtidas fazendo-se a média ao longo dos níveis do fator bloco
#          mu b1   b2   b3   b4   t1 t2 t3
rbind(t1=c(1, 1/4, 1/4, 1/4, 1/4, 1, 0, 0),
      t2=c(1, 1/4, 1/4, 1/4, 1/4, 0, 1, 0),
      t3=c(1, 1/4, 1/4, 1/4, 1/4, 0, 0, 1))%*%theta

# matrizes de incidência dos efeitos dos termos do modelo
X <- mapply(model.matrix, object=list(mu=~1, bl=~-1+bl, tr=~-1+tr),
            MoreArgs=list(data=da))

# as matrizes são irrestritas, ou seja, não são omitidas colunas
X <- do.call(cbind, X)

# valores esperados para cada observação do experimento
da$eta <- X%*%theta

# valor observado para cada observação do experimento
# uma variância pequena foi usada para compararmos os resultados
da$y <- da$eta+rnorm(nrow(da), 0, 0.01)

with(da, tapply(eta, list(bl, tr), function(x) x)) # esperados
with(da, tapply(y, list(bl, tr), function(x) x))   # observados

#-----------------------------------------------------------------------------
# estimação dos efeitos do modelo

m0 <- lm(y~bl+tr, data=da)
coef(m0) # estimativas dos efeitos, apenas k-1 de cada fator são estimados
with(da, tapply(y, tr, mean)) # médias amostrais marginais para tratamentos

# obtidas fazendo-se a média ao longo dos níveis do fator bloco
# estamos usando as estimativas e k-1 efeitos para cada fator
#          mu b2   b3   b4   t2 t3
rbind(t1=c(1, 1/4, 1/4, 1/4, 0, 0),
      t2=c(1, 1/4, 1/4, 1/4, 1, 0),
      t3=c(1, 1/4, 1/4, 1/4, 0, 1))%*%coef(m0)

require(gmodels)

# estimable() é uma forma diferente de fazer essa conta
#                         mu b2   b3   b4   t2 t3
estimable(m0, rbind(tr1=c(1, 1/4, 1/4, 1/4, 0, 0),
                    tr2=c(1, 1/4, 1/4, 1/4, 1, 0),
                    tr3=c(1, 1/4, 1/4, 1/4, 0, 1)))

# deve se ter habilidade para manipilar essas funções lineres
# elas mudam com o tipo de contraste empregado

require(contrast)
# contrast() é outra forma de fazer a conta que não exige conhecimento
# da parametrização usada, type="average" é para fazer média nos blocos
contrast(m0, type="average", list(bl=levels(da$bl), tr="t1"))
contrast(m0, type="average", list(bl=levels(da$bl), tr="t2"))
contrast(m0, type="average", list(bl=levels(da$bl), tr="t3"))

c0 <- contrast(m0, type="average", list(bl=levels(da$bl), tr="t3"))
c0$X # vetor de coeficientes, foi o que passamos para a estimable()

# independência entre efeitos de blocos e tratamentos
round(vcov(m0)*100,3)

#-----------------------------------------------------------------------------
# teste de Tukey para contraste entre médias (populacionais marginais)

require(agricolae)

# esse teste é baseado nas médias amostrais
with(da, HSD.test(y, tr, df.residual(m0), deviance(m0)/df.residual(m0)))

#-----------------------------------------------------------------------------
# o que fazer quando perdemos uma observação? entrar em pânico?
# todos os procedimento matriciais acima continuam válidos

da$z <- da$y
da$z[1] <- NA
with(da, tapply(z, list(bl, tr), function(x) x)) # 1 valor perdido

m1 <- lm(z~bl+tr, data=da)
coef(m1) # os efeitos ainda são estimáveis
round(vcov(m1)*100,6) # perde a ortogonalidade dos efeitos
with(da, tapply(z, tr, mean, na.rm=TRUE)) # médias amostrais são viesadas

# médias de mínimos quadrados são estimadores não viesados
# o desbalanceamento só causou diferença de precisão nas médias
#                         mu b2   b3   b4   t2 t3
estimable(m1, rbind(tr1=c(1, 1/4, 1/4, 1/4, 0, 0),
                    tr2=c(1, 1/4, 1/4, 1/4, 1, 0),
                    tr3=c(1, 1/4, 1/4, 1/4, 0, 1)))

contrast(m1, type="average", list(bl=levels(da$bl), tr="t1"))
contrast(m1, type="average", list(bl=levels(da$bl), tr="t2"))
contrast(m1, type="average", list(bl=levels(da$bl), tr="t3"))

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

# esse teste é baseado nas médias amostrais que agora são viesadas
# média harmônica o número de repetições e usado
with(da, HSD.test(z, tr, df.residual(m0), deviance(m0)/df.residual(m0)))

# viés das médias e aproximação do número de repetições -> teste não confiável

#-----------------------------------------------------------------------------
# contraste entre médias populacionais marginais

contrast(m1, type="average",
         list(bl=levels(da$bl), tr="t1"),
         list(bl=levels(da$bl), tr="t2"))
contrast(m1, type="average",
         list(bl=levels(da$bl), tr="t1"),
         list(bl=levels(da$bl), tr="t3"))
contrast(m1, type="average",
         list(bl=levels(da$bl), tr="t2"),
         list(bl=levels(da$bl), tr="t3"))

# essa sequência de testes de hipótese tem controle do nível nomínal de
# significância por teste, mas não controle global sobre todos os testes
# para isso deve se adotar o procedimento de Bonferroni ou outro
# não requer usar type="average" porque os efeitos de bloco se
# cancelam no contraste

#-----------------------------------------------------------------------------
# multcomp tem procedimentos que controlar o nível global de significância
# para o conjunto de hipóteses de interesse

require(multcomp)
summary(glht(m1, linfct=mcp(tr="Tukey")))

# Tukey não quer dizer teste de Tukey, e sim contrastes de Tukey, ou seja,
# todos as comparações duas a duas entre as médias

#-----------------------------------------------------------------------------
# níveis do fator bloco tem efeito aleatório

require(nlme)
m2 <- lme(z~tr, random=~1|bl, data=da, na.action=na.omit)

# não tem diferença de precisão nas estimativas das médias
estimable(m2, rbind(tr1=c(1,0,0), tr2=c(1,1,0), tr3=c(1,0,1)))
contrast(m2, list(tr=levels(da$tr)))

# contrastes devem ser feitos com as médias marginais sempre
# no caso de balanceamento, as médias amostrais são estimadores não viesados
# das médias populacionais marginais

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

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.

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