Arquivo

Posts Tagged ‘anova’

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()

#-----------------------------------------------------------------------------
Anúncios

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()

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

Método gráfico para valores iniciais em regressão não linear com RStudio

Obtenção de valores iniciais para o modelo duplo van Genuchten de 7 parâmetros com os recursos da função manipulate do RStudio.

No post Método gráfico interativo para valores iniciais em regressão não linear eu apresentei alguns procedimentos gráficos para obter bons valores iniciais para ajuste de modelos de regressão não linear. Ao obter valores iniciais pelo procedimento gráfico facilita-se a convergência do método de estimação e você passa a entender/interpretar melhor os parâmetros do modelo.

Nesse post vou apresentar os métodos gráficos construídos com os recursos do pacote manipulate que é parte do RStudio, um recente editor de scripts R que tem ganhado muitos adeptos por ser muito amigável. Eu já postei o uso dessas funções em RStudio e a função manipulate() onde falo sobre as opções das função.

Vou apresentar o procedimento com dados reais. Preciso encontrar chutes iniciais para o modelo duplo van Genuchten (ver artigo original). Esse modelo tem 7 parâmetros na seguinte equação

\theta(\Psi) = \theta_{res}+\displaystyle \frac{\theta_{pmp}-\theta_{res}}{(1+(\alpha_{tex}\Psi)^{n_{tex}})^{1-1/n_{tex}}}+\displaystyle \frac{\theta_{sat}-\theta_{pmp}}{(1+(\alpha_{est}\Psi)^{n_{est}})^{1-1/n_{est}}},

em que \theta e \Psi são valores observados e o restante são os parâmetros da função. Em geral, quanto maior o número de parâmetros a serem estimados num modelo de regressão mais demorara será a convergência porque o espaço da solução é de alta dimensão. Além do mais, com muitos parâmetros a serem estimados, pode haver fraca de identificabilidade do modelo. Dessa forma, bons chutes iniciais podem ajudar na convergência. Uma ligeira noção do campo de variação de cada parâmetro o usuário precisa ter (limites inferior e superior). Você pode inicialmente usar intervalos amplos e ir reduzindo o comprimento do intervalo a medida que observar quais são os valores mais próximos dos ótimos.

#-----------------------------------------------------------------------------
# dados para a Curva de Retenção de Água, umidade (theta) e tensão (psi)

cra <- data.frame(psi=c(0.01,1,2,4,6,8,10,33,60,100,500,1500,2787,4727,6840,
                    7863,9030,10000,10833,13070,17360,21960,26780,44860,
                    69873,74623,87287,104757,113817,147567,162723,245310,
                    262217,298223),
                  theta=c(0.779,0.554,0.468,0.406,0.373,0.36,0.344,0.309,
                    0.298,0.292,0.254,0.241,0.236,0.233,0.223,0.202,0.172,
                    0.187,0.138,0.098,0.07,0.058,0.052,0.036,0.029,0.0213,
                    0.0178,0.0174,0.0169,0.0137,0.0126,0.0109,0.0106,0.0053))

#-----------------------------------------------------------------------------
# gráfico dos valores observados

plot(theta~log10(psi), data=cra)

#-----------------------------------------------------------------------------
# função do modelo duplo van Genuchten, 7 parâmetros

dvg <- function(x, ts, ti, tr, ae, at, ne, nt){
  tr+(ti-tr)/((1+(at*10^x)^nt)^(1-1/nt))+(ts-ti)/((1+(ae*10^x)^ne)^(1-1/ne))
}

dvg(log10(c(1,10,100,1000,10000)),         # log 10 das tensões
    0.7, 0.2, 0.05, 1.3, 0.0001, 1.5, 3.5) # valor dos parâmetros

#-----------------------------------------------------------------------------
# procedimento gráfico para obter chutes e conhecer a função dos parâmetros

require(manipulate) # carrega pacote que permite manipulação gráfica
start <- list()     # cria uma lista vazia para receber os valores finais

manipulate({
             plot(theta~log10(psi), data=cra)
             curve(dvg(x, ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt), add=TRUE)
             start <<- list(ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt)
           },
           ts=slider(0.7, 0.9, initial=0.8),
           ti=slider(0.15, 0.25, initial=0.2),
           tr=slider(0, 0.10, initial=0.05),
           ae=slider(1.01, 3, initial=1.3),
           at=slider(0, 0.0001, initial=0.00005),
           ne=slider(1.01, 3, initial=1.65),
           nt=slider(1.8, 5, initial=4.3))

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

start # valores salvos do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo aos dados com os chutes do procedimento gráfico

n0 <- nls(theta~tr+(ti-tr)/((1+(at*psi)^nt)^(1-1/nt))+
          (ts-ti)/((1+(ae*psi)^ne)^(1-1/ne)),
          data=cra, start=start, trace=TRUE)
summary(n0) # quadro de estimativas

#-----------------------------------------------------------------------------
# gráfico dos dados com a curva estimada

lis <- c(list(x=NULL), as.list(coef(n0)), body(dvg))
plot(theta~log10(psi), cra,
     ylab=expression(Conteúdo~de~água~no~solo~(theta*","~g~g^{-1})),
     xlab=expression(Tensão~matricial~(Psi*","~kPa)), xaxt="n")
tmp <- as.function(lis)
curve(tmp, add=TRUE, col=2, lwd=1.5)
axis(1, at=-2:5, label=as.character(10^(-2:5)), lwd.ticks=2)
s <- log10(sort(sapply(1:9, function(x) x*10^(-3:6))))
axis(1, at=s, label=NA, lwd=0.5)
abline(v=-2:6, h=seq(0,1,0.05), lty=3, col="gray50")

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

O procedimento permite encontrar chutes iniciais para duas ou mais curvas. Basta escrever adequadamente as funções dentro da manipulate(). Abaixo estão os códigos para obter os valores inicias para o modelo logístico usado no contexto de epidemiologia. Em duas variedades de uma espécie vegetal foram observadas a severidade da doença ao longo do tempo. O objetivo é encontrar estimativas dos parâmetros para a função (que no caso é a logística) que associa severidade ao tempo e comparar as funções entre as variedades. A comparação das funções pode ser feita fazendo-se testes de hipóteses para funções dos parâmetros (em geral contrastes) ou testes envolvendo modelos aninhados, como você verá a seguir.

#-----------------------------------------------------------------------------
# dados de severidade em função do tempo para duas variedades

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

#-----------------------------------------------------------------------------
# gráficos exploratórios

require(lattice)
xyplot(V1+V2~DAI, data=da, type=c("p","a"))

# ajustar modelos. separado ou junto?
# qual modelo? quais os chutes?

#-----------------------------------------------------------------------------
# mudando a disposição dos dados

require(reshape)
db <- melt(da, id=c("DAI","rep"))
names(db)[3:4] <- c("vari","sever")
str(db)

#-----------------------------------------------------------------------------
# carrega pacote e faz a função para encontrar via deslizadores valores ótimos

require(manipulate)
start <- list()

manipulate({
  plot(sever~DAI, db, col=ifelse(db$vari=="V1",1,2))
  A1 <- a1; S1 <- s1; X01 <- x01
  A2 <- a2; S2 <- s2; X02 <- x02
  curve(A1/(1+exp(-(x-X01)/S1)), add=TRUE, col=1)
  curve(A2/(1+exp(-(x-X02)/S2)), add=TRUE, col=2)
  start <<- list(A=c(A1,A2), S=c(S1,S2), X0=c(X01,X02))
  },
  a1=slider(90, 110, initial=95, step=1),
  a2=slider(90, 110, initial=95, step=1),
  s1=slider(2, 8, initial=6, step=0.1),
  s2=slider(2, 8, initial=6, step=0.1),
  x01=slider(10, 40, initial=20, step=1),
  x02=slider(10, 40, initial=20, step=1)
  )

start # valores do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo não restrito: cada variedade com os seus parâmetros

n0 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S[vari])),
          data=db, start=start)
summary(n0)

#-----------------------------------------------------------------------------
# ajuste do modelo restrito: as variades possuem mesmo parâmetro S

start$S <- mean(start$S)
start

n1 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S)),
          data=db, start=start)
summary(n1)

#-----------------------------------------------------------------------------
# testa a hipótese do modelo nulo n1 tendo como alternativa o modelo n0

anova(n1, n0)

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

Obtenção de valores iniciais para o ajuste de duas curvas.

Ainda é necessário fazer a checagem das pressuposições envolvidas no modelo de regressão clássico: homocedasticidade (a variância dos desvios de regressão é constante ao longo a curva), normalidade (os erros aleatórios têm distribuição normal) e independência entre observações (normalmente a independência é garantida por um plano de amostragem adequado). Consulte Análise de resíduos em regressão não linear para saber como obter os 4 gráficos de análise de resíduos.

Esse post foi um aprimoramento de uma resposta na R-help e estimulado por citações na lista. Até a próxima ridícula.

Fatorial com desdobramento da interação, gráfico e tabela de médias

Rendimento de grãos em função dos níveis de potássio e umidade do solo. Letras minúsculas comparam níveis de potassio em um mesmo nível de umidade, as letras maiúsculas fazem o inverso.

O R não é um aplicativo específico de análise de experimentos. Mesmo assim, ele oferece funções que permitem a análise dos mais diversos tipos de experimentos, ficando a cargo do usuário saber usar corretamente as funções/pacotes que o R oferece.

Pode-se dizer que os experimentos fatoriais são os mais comuns na ciência. Nas ciências agrárias são amplamente usados. Têm-se uma tradição de representar a interação entre os fatores (qualitativos) empregando testes de médias para níveis de um fator fixando níveis dos demais.

Nessa ridícula, apresento procedimentos para obter o desdobramento de interação em um fatorial duplo. Os dados são do rendimento de grãos de soja em função dos fatores umidade do solo e dose de fertilizante potássico. O experimento foi instalado em casa de vegetação com as unidades experimentais (vasos) agrupadas em blocos que controlaram variações de luminosidade/umidade/ventilação.

Os procedimentos incluem o ajuste do modelo, o teste sobre os efeitos dos fatores pela análise de variância, a checagem das pressuposições do modelo, o fatiamento das somas de quadrados, o teste de médias. No final são apresentados procedimentos para representar os resultados em gráfico e tabela de dupla entrada. É importante dizer que os fatores deveriam ser analisados como quantitativos, onde espera-se uma interpretação mais interessante sobre o fenômeno. Entretanto, o método abaixo é ilustrativo do procedimento para fatores qualitativos. Em uma futura oportunidade analisarei dos dados considerando os fatores como quantitativos.

#-----------------------------------------------------------------------------
# importa e prepara os dados
 
rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt",
                   header=TRUE, sep="\t")
rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
str(rend)
 
#-----------------------------------------------------------------------------
# análise gráfica
 
require(lattice)
xyplot(rg~K|A, groups=bloc, data=rend, type="b", auto.key=TRUE)
 
#-----------------------------------------------------------------------------
# ajuste do modelo e faz checagem
 
m0 <- aov(rg~bloc+A*K, data=rend)
summary(m0)
par(mfrow=c(2,2)); plot(m0); layout(1)
 
#-----------------------------------------------------------------------------
# desdobrando somas de quadrados para a variação de K dentro de A e vice-versa

m1 <- aov(rg~bloc+K/A, data=rend)
desAinK <- sapply(paste("K", levels(rend$K), sep=""), simplify=FALSE,
                  grep, x=names(coef(m1)[m1$assign==3]))
summary(m1, split=list("K:A"=desAinK))

m2 <- aov(rg~bloc+A/K, data=rend)
desKinA <- sapply(paste("A", levels(rend$A), sep=""), simplify=FALSE,
                  grep, x=names(coef(m2)[m2$assign==3]))
summary(m2, split=list("A:K"=desKinA))
 
#-----------------------------------------------------------------------------
# desdobrando a interação em testes de Tukey
 
require(agricolae)
require(plyr)

KinA <- sapply(levels(rend$A), simplify=FALSE,
               function(a){
                 with(subset(rend, A==a),
                      HSD.test(rg, K,
                               DFerror=df.residual(m0),
                               MSerror=deviance(m0)/df.residual(m0)))
               })

KinA <- ldply(KinA, NULL)
KinA$M <- gsub(" ", "", KinA$M, fixed=TRUE)
KinA$trt <- as.factor(as.numeric(as.character(KinA$trt)))
str(KinA)

AinK <- sapply(levels(rend$K), simplify=FALSE,
               function(k){
                 with(subset(rend, K==k),
                      HSD.test(rg, A,
                               DFerror=df.residual(m0),
                               MSerror=deviance(m0)/df.residual(m0)))
               })

AinK <- ldply(AinK, NULL)
AinK$M <- toupper(gsub(" ", "", AinK$M, fixed=TRUE))
AinK$trt <- as.factor(as.numeric(as.character(AinK$trt)))
str(AinK)

#-----------------------------------------------------------------------------
# combinando os data.frames para obter o gráfico e a tabela

KinA$comb <- paste(KinA$trt, KinA$.id)
AinK$comb <- paste(AinK$.id, AinK$trt)
desd <- merge(KinA, AinK, by=c("comb","comb"))
desd$let <- with(desd, paste(M.y, M.x, sep=""))
desd <- desd[,c("trt.x","trt.y","means.x","let")]
names(desd) <- c("K","A","means","let")
str(desd)

#-----------------------------------------------------------------------------
# gráfico das médias com as letras sobre as barras

#png("f005.png", w=500, h=250)
barchart(means~K|A, data=desd, horiz=FALSE, layout=c(3,1),
         ylab="Rendimento de grãos (g/vaso)", ylim=c(NA,45),
         xlab="Dose de potássio (mg/dm³)",
         key=list(text=list("Níveis de umidade do solo (%)")),
         sub=list("*letras minúsculas comparam médias em um mesmo nível de umidade.", cex=0.8, font=3),
         strip=strip.custom(bg="gray90"), col=colors()[51],
         panel=function(x, y, subscripts, ...){
           panel.barchart(x, y, subscripts=subscripts, ...)
           panel.text(x, y, label=desd[subscripts,"let"], pos=3)
         })
#dev.off()

#-----------------------------------------------------------------------------
# tabela de dupla entrada com as médias das combinações A e K

require(reshape)
desd$means.let <- with(desd, paste(means, let))
cast(desd, K~A, value="means.let")

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

Ajuste de muitos modelos de regressão linear

Pessoas que estão migrando de outro aplicativo para o R normalmente vivem perguntando “o R faz isso?, o R faz aquilo?…”. A resposta é “Sim, o R faz!”. Circula entre usuários o seguinte ditado: “não se pergunta o que o R faz, e sim como ele faz”.

Acontece que o R não é um software específico de ajuste de modelos de regressão, análise multivariada, análise de experimentos ou qualquer outra técnica. O R é geral. Os pacotes do R, por sua vez, são capazes de concentrar funções específicas. Um aplicativo que só faz ajuste de modelos de regressão tem espaço de sobra para incluir diversas opções, contem lista modelos para ajustar aos dados, enfim, mas tudo isso porque é destinado apenas para esse fim.

Por outro lado, no R é possível desenvolver procedimentos para qualquer tipo de tarefa. Se você quer ajustar uma lista de modelos e depois escolher o de maior R², ótimo! No você também consegue isso. Veja no CMR abaixo.

# gera dados
da <- data.frame(x=runif(100), z=5*rpois(100, lambda=7), w=runif(100, 50, 100))
da$y <- with(da, 12+0.1*x+0.05*z+0.34*w+0.2*sqrt(z)+0.1*x*w)+rnorm(100,0,0.1)
 
# vetor com as fórmulas específicando diferentes modelos lineares
form <- c(mod1=y~x, mod2=y~x+z, mod3=y~x+I(x^2), mod4=y~x+z+w)
 
# ajuste dos modelos
ajustes <- lapply(form, function(f){ m0 <- lm(f, data=da); m0 })
 
lapply(ajustes, summary) # quadro geral de estimativas e qualidade
lapply(ajustes, anova)   # quadro de anova sequencial
lapply(ajustes, coef)    # vetor de estimativas
sapply(ajustes, function(a){ summary(a)$r.squared})     # R²
sapply(ajustes, function(a){ summary(a)$adj.r.squared}) # R² ajustado
sapply(ajustes, function(a){ summary(a)$sigma})         # QMR
sapply(ajustes, deviance)                               # SQR
sapply(ajustes, df.residual)                            # GLR
lapply(ajustes, function(a){ summary(a)$coeff})         # tabela de estimativas
do.call(rbind, lapply(ajustes, function(a){ summary(a)$coeff})) # junta das tabelas
sapply(ajustes, fitted)    # valores ajustados
sapply(ajustes, residuals) # resíduos da análise
sapply(ajustes, vcov)      # matriz de covariância das estimativas
apply(sapply(ajustes, residuals), 2, shapiro.test) # normalidade dos resíduos

Sempre faça a avaliação das pressuposições do modelo antes de aplicar inferências. A escolha de modelo apenas pelo valor do R² não é aconselhada.

Realçador de código R para blogs

Quem trabalha com programação reconhece que a organização do código facilita muito a releitura e checagem. O R, uma vez que é uma linguagem, possui editores que apresentam recursos para organização do código, como Emacs, RStudio, Tinn-R, Rkward entre diversos. A maioria dos editores apresenta um padrão de identação do código e realçador (highlighter) de texto. A necessidade de divulgação de código R, de forma organizada e atraente, fez com que os realçadores aparecessem nas páginas de blogs, wikis, etc. O wordpress.com possui um realçador de código R nativo. Para quem não tem wordpress.com, pode usar o Pretty R syntax highlighter que possui, além dos realces usuais, links para as páginas de documentação das funções. Abaixo está código R com o realçador do wordpress.com e com o pretty-R. Perceba que as funções (negrito azul) possuem links para página de ajuda.

# importa dados
soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)
 
# ajusta um modelo e pede anova
m1 <- aov(rengrao~bloco+agua*potassio, soja)
anova(m1)
 
# cria uma lista com as variáveis resposta
respostas <- do.call(c, apply(soja[,4:7], 2, list))
do.call(c, respostas)
 
# faz o ajuste para todas as respostas
ajustes <- lapply(respostas,
                  function(r){
                    m0 <- aov(r~bloco+agua*potassio, data=soja)
                    m0
                  })
 
# pede todas as anovas
lapply(ajustes, anova)
# importa dados
soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)
 
# ajusta um modelo e pede anova
m1 <- aov(rengrao~bloco+agua*potassio, soja)
anova(m1)
 
# cria uma lista com as variáveis resposta
respostas <- do.call(c, apply(soja[,4:7], 2, list))
do.call(c, respostas)
 
# faz o ajuste para todas as respostas
ajustes <- lapply(respostas,
                  function(r){
                    m0 <- aov(r~bloco+agua*potassio, data=soja)
                    m0
                  })
 
# pede todas as anovas
lapply(ajustes, anova)

Created by Pretty R at inside-R.org

Além de no R você fazer praticamente tudo, você pode apresentar de forma elegante. Até a próxima ridícula.