Arquivo

Posts Tagged ‘lattice’

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

#-----------------------------------------------------------------------------
Anúncios
Categorias:delineamento Tags:, ,

Legenda e contornos em gráficos de nível

Peso de 100 grãos em função do nível de saturação de água e potássio. Contornos e rótulo na escala de cores são os destaques desse gráfico.

Mais uma ridícula de lattice. Fiz duas coisas: adicionei rótulo à legenda de cores (colorkey) e coloquei contornos sobre o gráfico de níveis. Assim como a maioria das dicas gráficas, essa também é baseada nas mensagens da r-help. Os respectivos links estão no CMR. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# leitura dos dados

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

#-----------------------------------------------------------------------------
# ver o peso de 100 grãos

require(lattice)

xyplot(pesograo~potassio, groups=agua, data=soja, type=c("p","a"))
xyplot(pesograo~potassio|agua, data=soja, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste de um modelo polinomial nos dois fatores

m0 <- lm(pesograo~bloco+poly(agua,2)*poly(potassio,2), data=soja)
par(mfrow=c(2,2)); plot(m0); layout(1)

summary(influence.measures(m0))
soja[55,] # influente segundo dffits

m0 <- lm(pesograo~bloco+poly(agua,2)*poly(potassio,2),
         data=soja[-55,])
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
summary(m0)

m1 <- lm(pesograo~bloco+(agua+I(agua^2))*potassio+I(potassio^2),
         data=soja[-55,])
par(mfrow=c(2,2)); plot(m1); layout(1)

anova(m0, m1)
summary(m1)

#-----------------------------------------------------------------------------
# fazendo a predição

pred <- expand.grid(bloco="I", agua=seq(37.5,62.5,l=30),
                    potassio=seq(0,180,l=30))
pred$y <- predict(m1, newdata=pred)

#-----------------------------------------------------------------------------
# representando com wireframe()

require(RColorBrewer)

display.brewer.all()
colr <- brewer.pal(11, "RdYlGn")
colr <- colorRampPalette(colr, space="rgb")

zlab <- "Peso de 100 grãos"
xlab <- "Potássio no solo"
ylab <- "Nível de saturação de água"

wireframe(y~potassio*agua, data=pred,
          scales=list(arrows=FALSE), zlab=list(zlab, rot=90),
          xlab=list(xlab, rot=24), ylab=list(ylab, rot=-37),
          col.regions=colr(100),  drape=TRUE,
          screen=list(z=40, x=-70))

#-----------------------------------------------------------------------------
# representando em um levelplot()

# grid mais fino
pred <- expand.grid(bloco="I", agua=seq(37.5,62.5,l=100),
                    potassio=seq(0,180,l=100))
pred$y <- predict(m1, newdata=pred)

levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
          xlab=xlab, ylab=ylab)

#-----------------------------------------------------------------------------
# adicionando rotulo à legenda de cores, baseado nas mensagens da r-help
# http://r.789695.n4.nabble.com/Adding-title-to-colorkey-td4633584.html

library(grid)

# modo 1
levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
          xlab=xlab, ylab=ylab,
          par.settings=list(layout.widths=list(axis.key.padding=4)))
grid.text(zlab, x=unit(0.88, "npc"), y=unit(0.5, "npc"), rot=90)

# modo 2
levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
          xlab=xlab, ylab=ylab,
          ylab.right=zlab,
          par.settings=list(
            layout.widths=list(axis.key.padding=0, ylab.right=2)))

require(latticeExtra)

# modo 3
p <- levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
               xlab=xlab, ylab=ylab,
               par.settings=list(
                 layout.widths=list(right.padding=4)))
p$legend$right <- list(fun=mergedTrellisLegendGrob(p$legend$right,
                         list(fun=textGrob, args=list(zlab, rot=-90, x=2)),
                         vertical=FALSE))
print(p)

#-----------------------------------------------------------------------------
# adicionando contornos, baseado em
# https://stat.ethz.ch/pipermail/r-help/2006-February/088166.html

#png(file="f037.png", width=500, height=400)
p <- levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
               xlab=xlab, ylab=ylab,
               panel=function(..., at, contour=FALSE, labels=NULL){
                 panel.levelplot(..., at=at, contour=contour, labels=labels)
                 panel.contourplot(..., at=at, contour=TRUE,
                                   labels=list(labels=format(at, digits=4),
                                     cex=0.9))
               },
               par.settings=list(
                 layout.widths=list(right.padding=4)))
p$legend$right <- list(fun=mergedTrellisLegendGrob(p$legend$right,
                         list(fun=textGrob, args=list(zlab, rot=-90, x=2)),
                         vertical=FALSE))
print(p)
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:,

Gráfico de densidade com destaque dos quantis

Gráficos de densidade com regiões destacando os quartis.

Gráficos de densidade são utilizados para se conhecer a forma da distribuição dos dados. A função lattice::densityplot() faz esse tipo de gráfico. Para aprimorar a visualização dos quartis dos dados, eu implementei uma função (my.densityplot()) para que destaque com regiões coloridas as porções de dados delimitadas pelos quartis. O código abaixo gera a figura no início desse post. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# densityplot() com regiões pintadas delimitadas pelos quantis

require(lattice)

my.densityplot <- function(x, cols="gray90", probs=c(1,3)/4, ...){
  panel.densityplot(x, ...)
  dx <- density(x)
  breaks <- sort(c(range(dx$x), quantile(x, probs=probs)))
  fx <- approxfun(dx$x, dx$y)
  do.polygon <- function(x){
    y <- fx(x)
    return(list(x=c(min(x), x, max(x)), y=c(0, y, 0)))
  }
  seqs <- lapply(1:(length(breaks)-1),
                 function(i){
                   x <- seq(breaks[i], breaks[i+1], l=20)
                   do.polygon(x)
                 })
  for(i in 1:length(seqs)){
    seqs[[i]]$col <- cols[i]
  }
  lapply(seqs, function(i) do.call(panel.polygon, i))
  panel.mathdensity(dmath=dnorm, col="black",
                    args=list(mean=mean(x),sd=sd(x)), lty=3)
}

require(RColorBrewer)
display.brewer.all()

cols <- brewer.pal(3, "BuPu") # Greens, Blues, BuPu, Reds
densityplot(~height|voice.part, data=singer,
            cols=cols, plot.points="rug", panel=my.densityplot)

#png("f036.png", width=500, height=320)
cols <- brewer.pal(4, "Greens") # Greens, Blues, BuPu, Reds
densityplot(~height|voice.part, data=singer, layout=c(4,2),
            strip=strip.custom(bg="gray90"),
            xlab="Altura", ylab="Densidade", col=1,
            cols=cols, probs=1:3/4, plot.points="rug", panel=my.densityplot)
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:,

Equação junto ao gráfico na lattice

Reposta à salinidade para quatro espécies de craca. A regulação osmótica está ligada à capacidade adaptativa do organismo ao ambiente. O gráfico mostra os valores observados e a curva de regressão estimada acompanhada de bandas de confiança (95%). As equações em cada gráfico representam o modelo ajustado que descreve a curva.

Os polinômios são os modelos de regressão mais empregados em algumas áreas do conhecimento. Isso não quer dizer que eles sejam modelos bons. Na verdade, isso depende do cada um define como um bom modelo. Os polinômios são uma aproximação local. Isso significa que eles não devem ser usados para predições fora da amplitude observada das covariáveis (região experimental). Em outras palavras, não se deve fazer extrapolações, mas pode-se, com certos cuidados, fazer interpolações. Os polinômios são desprovidos de interpretação a partir do segundo grau. A partir do segundo grau, apenas o intercepto tem interpretação pois é o valor da função quando a covariável (x) é zero. Por isso, interpretações em geral são perdas de tempo. O que pode-se fazer são julgamentos quando ao sinal do coeficiente de maior grau que por exemplo, num modelo quadrático, indica se a curvatura é para cima (ou não). Outro ponto é que os polinômios apresentam problema de colinearidade. Para superar essa deficiência pode-se centrar a covariável ou usar polinômios ortogonais. Apesar do polinômio já ser sem interpretação, o uso de polinômios ortogonais é menos interpretável ainda pois os coeficientes não são taxas em relação à x, x², x³…, mas para uma transformação linear de x, essa que induz a ortogonalidade. Resumindo, os polinômios ortogonais são úteis pois minimizam o problema de colineridade. Então, do ponto de vista de reportar uma equação de predição, aquela que descreve o valor esperado de y (resposta) em função de x, é mais útil usar o polinômio cru. Isso porque, mesmo sem ter interpretação, toda pessoa sabe traçar a curva de um polinômio, é algo familiar, visto desde o ensino médio, e novamente, não que isso seja um ponto positivo para um modelo. De qualquer forma, o post de hoje ensina a ajustar polinômios e representar no gráfico os valores observados, os valores ajustados com intervalo de confiança e a equação dentro do gráfico. Os procedimentos abaixo eu desenvolvi para atender demandas de consultoria. Até a próxima ridícula.

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

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

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

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

xyplot(Osm~Salinidade|Espécie, data=osmo, type=c("p","smooth"))

#-----------------------------------------------------------------------------
# vou desconsiderar o fato dos dados serem heterocedasticos e ajustar um lm
# em outra matéria trataremos os desvios de pressupostos
# ajuste usando polinômios ortogonais

m0 <- aov(Osm~Espécie*poly(Salinidade,2), data=osmo)
par(mfrow=c(2,2)); plot(m0); layout(1)

qqmath(~residuals(m0)|osmo$Espécie)
xyplot(sqrt(abs(residuals(m0)))~fitted(m0)|osmo$Espécie, type=c("p","smooth"))
# desvio de pressupostos, trataremos em outra matéria do blog

summary(m0) # inferências sem validade pois existe desvio de pressupostos

#-----------------------------------------------------------------------------
# partição das somas de quadrados, declara Salinidade dentro de Espécie

m1 <- aov(Osm~Espécie/poly(Salinidade,2), data=osmo)

names(coef(m1))
attr(terms(m1), "term.labels")
ncoef <- 2 # número de efeitos por expécie=2 (linear e quadrático)
ngrup <- 4 # número de níveis de espécie=4
split.list <- split(matrix(1:(ncoef*ngrup), ngrup, ncoef), f=1:ngrup)
names(split.list) <- levels(osmo$Espécie)
split.list <- lapply(split.list,
                     function(x){ names(x) <- c("lin","qua"); x})
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

split.list <- as.list(1:(ncoef*ngrup))
names(split.list) <- c(outer(levels(osmo$Espécie), c("lin","qua"),
                             paste, sep=": "))
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

#-----------------------------------------------------------------------------
# prepara as equações para colocar no gráfico, função para polinômiais

coefs <- coef(lm(Osm~-1+Espécie/(Salinidade+I(Salinidade^2)), data=osmo))
coefs <- as.numeric(formatC(coefs, format="f", digits=3, flag="-"))

get.expression <- function(coef, digits){
  coef <- as.numeric(formatC(coef, format="f", digits=digits, flag="-"))
  postext <- c("","*x", paste("*x^", 2:(length(coef)-1), sep=""))
  sig <- c(ifelse(coef[1]<0, "-", ""), ifelse(coef[-1]<0, "-", "+"))
  txt <- paste(paste(sig, abs(coef), postext, sep=""), collapse=" ")
  eval(parse(text=paste("expression(y==", txt, ")", sep="")))
}

coefs <- split(matrix(coefs, ncol=3), f=1:4)
eqn <- lapply(coefs, get.expression, digits=3)
eqn

#-----------------------------------------------------------------------------
# nomes científicos em itálico nas tarjas

fl <- as.expression(lapply(levels(osmo$Espécie),
    function(x){ bquote(italic(.(x))) }))
fl

#-----------------------------------------------------------------------------
# predição

pred <- expand.grid(Espécie=levels(osmo$Espécie), Salinidade=seq(0,40,by=2))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))

#-----------------------------------------------------------------------------
# carrega funções que serão necessárias para fazer as bandas de confiança

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
grep("ciH", ls(), value=TRUE) # nome das funções (painéis) que vamos usar

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

#png("f035.png", width=600, height=500)
xyplot(Osm~Salinidade|Espécie, data=osmo, col=1,
       xlab=expression(Salinidade~(ppm)),
       ylab=expression(Osmolalidade~(mOsm~kg^{-1}~de~H[2]*O)),
       strip=strip.custom(bg="gray90", factor.levels=fl))+
  as.layer(xyplot(fit~Salinidade|Espécie, data=pred, type="l", col=1,
                  ly=pred$lwr, uy=pred$upr,
                  prepanel=prepanel.ciH, panel=panel.ciH))
dimLay <- trellis.currentLayout()
for(col in 1:ncol(dimLay)){
  for(row in 1:nrow(dimLay)){
    trellis.focus("panel", col, row, highlight=FALSE)
    grid.text(label=eqn[[dimLay[row,col]]],
              x=0.03, y=1-0.03, just=c("left","top"))
  }
}
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, ,

Gráfico tridimensional com contornos de nível

Peso de 100 grãos em função do nível de saturação de água no solo e conteúdo de potássio fornecido pela adubação. Linhas no plano inferior representam a projeção das linhas de mesmo valor para peso de 100 grãos.

Nesse post vou apresentar uma função para adicionar contornos de nível em gráficos tridimensionais gerados pela lattice::wireframe(). Para ilustrar o uso da função eu fiz um ajuste de modelo de regressão múltipla aos dados do experimento com soja. Esses dados eu já usei em várias matérias aqui no Ridículas. Uma vez ajustado o modelo, os valores preditos são então representados pela wireframe() usando o panel.3d.contour que foi uma função que aprimorei das leituras que fiz das mensagens da r-help.

Um gráfico tridimensional como este impressiona quem o vê, porém esteja ciente das más impressões que um gráfico tridimensional pode dar pois ele não é de fato tridimensional e sim uma projeção tridimensional em um plano, o que acarreta algumas deformações. Todavia, achei por bem documentar e expor a função mas isso não significa que esta é a melhor representação. A função permite colocar os contornos no plano abaixo, acima e sobre a superfície (bottom, top, on). Tenha cuidado ao escolher o ângulo de observação do gráfico (screen=) e controle a amplitude do eixo z (zlim=) para que a superfície não esconda a projeção abaixo. Confira os resultados no CMR abaixo. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# leitura dos dados

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

#-----------------------------------------------------------------------------
# ver o peso de 100 grãos

require(lattice)

xyplot(pesograo~potassio, groups=agua, data=soja, type=c("p","a"))
xyplot(pesograo~potassio|agua, data=soja, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste de um modelo polinomial nos dois fatores

m0 <- lm(pesograo~bloco+poly(agua,2)*poly(potassio,2), data=soja)
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
summary(m0)

m1 <- lm(pesograo~bloco+(agua+I(agua^2))*potassio+I(potassio^2), data=soja)
par(mfrow=c(2,2)); plot(m1); layout(1)

anova(m0, m1)
summary(m1)

#-----------------------------------------------------------------------------
# fazendo a predição

pred <- expand.grid(bloco="I", agua=seq(37.5,62.5,l=30),
                    potassio=seq(0,180,l=30))
pred$y <- predict(m1, newdata=pred)

#-----------------------------------------------------------------------------
# vendo a superfície ajustada

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
args(panel.3d.contour)
body(panel.3d.contour)

require(RColorBrewer)

colr <- brewer.pal(11, "Spectral")
colr <- colorRampPalette(colr, space="rgb")

zlab <- "Peso de 100 grãos"
xlab <- "Potássio no solo"
ylab <- "Nível de saturação de água"

#png(file="f034.png", width=500, height=500)
wireframe(y~potassio*agua, data=pred,
          scales=list(arrows=FALSE), zlab=list(zlab, rot=90),
          xlab=list(xlab, rot=24), ylab=list(ylab, rot=-37),
          zlim=c(5,16), col="gray50", col.contour=1,
          panel.3d.wireframe="panel.3d.contour", type=c("on","bottom"),
          col.regions=colr(100),  drape=TRUE,
          screen=list(z=40, x=-70))
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:,

Média e desvio-padrão de muitas variáveis separado por grupos

Gráfico de barras com as médias e desvios-padrões para variáveis resposta de um experimento avaliando cultivares de banana e meios de cultura em cultura de tecidos.

Em estatística descritiva é comum a síntese dos dados em medidas de posição e dispersão. As medidas mais frequentemente usadas são a média e o desvio-padrão. Não raramente, existe a necessidade de obter essas medidas separando as observações por alguma variável categórica do conjunto de dados. Nesse post apresento algumas formas de obter isso. Estou certo de que as soluções que apresento abaixo não cobrem 5% das disponíveis considerando a extensa quantidade de pacotes/funções do R ligadas à análise descritiva, portanto, considere isso apenas como um ponto de partida. Por fim, sintetizo em um gráfico os resultamos aplicados em um conjunto de dados reais. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# gerando um conjunto de dados artificial

da <- expand.grid(A=1:3, B=1:2, r=1:15)
da <- cbind(da, as.data.frame(matrix(rnorm(nrow(da)*3), nrow(da))))
str(da)

#-----------------------------------------------------------------------------
# usando funções do reshape e plyr

require(reshape)
require(plyr)

ddply(da[,-c(1:3)], .(da$A, da$B), mean) # as médias
ddply(da[,-c(1:3)], .(da$A, da$B), sd)   # os desvios-padrões

db <- melt(da, id.vars=1:3, id.measures=4:6)
str(db)

#-----------------------------------------------------------------------------
# fazendo empilhamento das observações

ddply(db[,-3], .(A, B, variable), summarise, m=mean(value), sd=sd(value))

#-----------------------------------------------------------------------------
# usando funções do pacote doBy

library(doBy)

summaryBy(V1+V2+V3~A+B, data=da,
          FUN=function(x){ c(m=mean(x), sd=sd(x)) })

#-----------------------------------------------------------------------------
# usando funções do Hmisc

require(Hmisc)

with(da, bystats(cbind(V1, V2, V3), A, B,
                 fun=function(x){ c(m=colMeans(x), s=sd(x)) }))

#-----------------------------------------------------------------------------
# dados reais

bnn <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/banana.txt",
                  header=TRUE, sep="\t")
str(bnn)
bnn <- melt(bnn, id.measures=4:12)
bnn <- ddply(bnn, .(cult, meio, variable),
             summarise, m=mean(value), sd=sd(value))
str(bnn)

#-----------------------------------------------------------------------------
# gráfico de barras com desvio-padrão

require(lattice)
require(RColorBrewer)
display.brewer.all()
cols <- brewer.pal(nlevels(bnn$meio), "Set3")

png("f028.png", 500, 500)
trellis.par.set(superpose.polygon=list(col=cols),
                strip.background=list(col="gray90"))
barchart(m~cult|variable, groups=meio, data=bnn,
         sdv=bnn$sd, scales="free", layout=c(3,3),
         xlab="Cultivar", ylab=expression(bar(x)+sd),
         key=list(title="Meio de cultura", cex.title=1.2,
           columns=3, text=list(levels(bnn$meio)),
           rectangle=list(col=trellis.par.get("superpose.polygon")$col)),
         prepanel=function(y, sdv, subscripts=subscripts, ...){
           ly <- as.numeric(y+sdv[subscripts])
           list(ylim=range(y, ly, finite=TRUE))
         },
         panel=function(x, y, subscripts, groups, sdv, box.ratio, ...){
           panel.barchart(x, y, subscripts=subscripts,
                          groups=groups, box.ratio=box.ratio, ...)
           d <- 1/(nlevels(groups)+nlevels(groups)/box.ratio)
           g <- (as.numeric(groups[subscripts])-1); g <- (g-median(g))*d
           panel.arrows(as.numeric(x)+g, y, as.numeric(x)+g, y+sdv[subscripts],
                        angle=90, length=0.025)
         })
dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, , ,

Uma alternativa para apresentar bandas de confiança em regressão

Valores preditos para índice agronômico de milho em função da dose de nitrogênio e cultivar. Regiões em cinza são os intervalos de confiança (95%) para os valores preditos.

Neste post novamente vou apresentar a análise dos dados de índice agronômico de milho que iniciei em outros posts (Predição em modelos de regressão com blocos, Como fazer regressão polinomial com diferentes graus). A única novidade será a representação gráfica, que ao invés de separar as curvas por nível de cultivar, vou apresentar no mesmo gráfico usando uma região escura para representar os intervalos de confiança para os valores preditos. Esse gráfico facilita perceber se os intervalos têm sobreposição ou não. No entanto, ter (ou não) sobreposição não aponta aceitação (ou rejetição) da hipótese de igualdade dos valores preditos. Hipóteses devem ser definidas antes da realização do experimento e avaliadas por testes adequados. Foi necessário desenvolver uma função complemento para os gráficos da lattice para produzir as regiões de confiança. No script abaixo eu reproduzo as análises e faço o gráfico que mencionei. Até a próxima ridícula!

#-----------------------------------------------------------------------------
# dados de índice agronômico de milho

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

#-----------------------------------------------------------------------------
# ajuste do modelo. veja os posts anteriores a esse mencionados no texto.

levels(da$cultivar)
da$Pi.C <- with(da, ifelse(da$cultivar=="Pioneer-B815", dose^3, 0))
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)

# vetor de estimativas sem os efeitos de blocos
b <- coef(m3)[-grep("bloco", names(coef(m3)))]

# margem de erro, produto do erro padrão e quantil
me <- predict(m3, newdata=cbind(bloco=levels(da$bloco)[1], pred),
              se.fit=TRUE)$se*qt(0.975, df=df.residual(m3))

# data.frame com valores preditos e IC
pred$fit <- X%*%b
pred$lwr <- pred$fit-me; pred$upr <- pred$fit+me
str(pred)

#-----------------------------------------------------------------------------
# paineis feitos para lattice para apresentar os intervalos de confiança como
# regiões de outra cor, úteis para sobrepor ICs

prepanel.ciH <- function(x, y, ly, uy, subscripts, ...){
  x <- as.numeric(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  list(ylim=range(uy, ly, finite=TRUE), xlim=range(x))
}

panel.ciH <- function(x, y, ly, uy, subscripts, ...){
  y <- as.numeric(y)
  x <- as.numeric(x)
  or <- order(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  panel.polygon(c(x[or], x[rev(or)]),
                c(ly[or], uy[rev(or)]),
                col=1, alpha=0.03, border=NA)
  panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
  panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
  panel.xyplot(x, y, subscripts=subscripts, ...)
}

#-----------------------------------------------------------------------------
# gráfico com observados, preditos e IC pelo método matricial

require(lattice)
require(latticeExtra) # as.layer()
require(grid)

#png("f020.png", w=500, h=350)
with(pred,
     xyplot(fit~dose, groups=cultivar, type="l", lwd=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),
            ly=lwr, uy=upr,
            prepanel=prepanel.ciH,
            panel=panel.superpose,
            panel.groups=panel.ciH,
            key=list(space="top", columns=1,
              points=list(pch=1), text=list("valores observados"),
              lines=list(col=1), text=list("valores preditos"),
              rectangle=list(col="gray90", border=NA), text=list("IC 95%"))
            ))+
  as.layer(xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE))
draw.key(list(columns=1,
              lines=list(lty=1, lwd=2,
                col=trellis.par.get("superpose.symbol")$col[1:3]),
              text=list(levels(da$cultivar))),
         draw=TRUE, vp=viewport(x=unit(0.75, "npc"), y=unit(0.25, "npc")))
#dev.off()

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