Arquivo

Posts Tagged ‘fatorial’

Fatorial duplo com dois tratamentos adicionais

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

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

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

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

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

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

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

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

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

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

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

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

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

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

xtabs(~fungicida+aplicacao, fu)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

str(cima)
cima$confint

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Fuzzy c-means, um exemplo!

Pois bem, por sugestão do Professor Julio da Motta Singer apresentarei um script para realizar uma análise de agrupamento no escopo da lógica fuzzy. Seja um experimento fatorial (a mesma natureza dos dados aqui disponibilizados) a seguinte implementação realiza a clusterização de uma dado fator marginalmente ao outro fator.

A heurística da lógica fuzzy tem um incremento interpretativo de, ao contrário do agrupamento k-means, que é dicotômico, associar a cada cluster uma “probabilidade” de associação. Todavia, no método arbitramos uma constante (m) de fuzzyficação. Assim o seguinte script realiza o agrupamento fuzzy (por intermédio da função cmeans() da biblioteca e1071) para vários valores m assim é possível observar o comportamento (decaimento) do erro associado a clusterização e também a disposição dos fatores em relação ao cluster.

Infelizmente (ou para alegria geral) não vou me atrever a discorrer da teoria envolvida nessa análise e um estudo específico deixo a cargo dos maiores interessados.

Todo esse mise en cene, me soa um pouco genérico demais, eu sei (–Me desculpem!). Mas acredito que o script possa ser útil de alguma forma.

Aproveito o ensejo para agradecer as sugestões do Prof Julio da Motta Singer ao trabalho e o incentivo a esse post. Também deixo um alô ao colegas que participaram, de uma forma ou de outra!

Vamos ao trabalho, primeiro instalando (para quem ainda não o fez) e carregando a biblioteca e o conjunto de dados!

## install.packages('e1071')
library(e1071) # carrega pacote com implementação 'c-means'

dados <- transform(read.table("http://dl.dropbox.com/u/38195533/fac_exemplo.txt", # arquivo
                              header = TRUE,
                              sep = '\t',
                              na.string = '.',
                              dec = ','),
                   fatorA = factor(fatorA),
                   fatorB = factor(fatorB))

str(dados)

ajuste <- lm(resp ~ fatorA + fatorB, data = dados)
anova(ajuste) # ANOVA

matriz <- tapply(dados$resp, list(dados$fatorB, dados$fatorA), mean)

A primeira etapa foi realizada. Não posso passar ou próximo trecho sem deixa-los a par do quem está por vir! O agrupamento que “tentamos” é com relação ao um padrão de comportamento do fator sob análise, em relação ao outro fator, assim constituímos ideótipos e assim procedemos. O modo como constituimos esses ideótipos é aplicado unicamente a natureza desse conjunto de dados, mas outros ideótipos, ao seu critério podem ser usados, fique a vontade!

Antes de passarmos agora a análise em si manipulamos os dados e constituímos os centróides:

FinWil <- apply(matriz, 2, mean) - mean(matriz) # Finlay & Wilkinson

altos <- which(FinWil >= 0) # ambientes acima da media
baixo <- which(FinWil < 0) # ambientes abaixo da media

medias <- cbind(apply(matriz[, altos], 1, mean),
                apply(matriz[, baixo], 1, mean))

##--------------------------------------------------------------------
## IDEOTIPOS
## I   - maximo em todos 
## II  - maximo nos altos minimo nos baixos
## III - minimo nos altos maximo nos baixos
## IV  - minimo em todos 
## ADICIONAIS:
## V   - media em todos
## VI  - maximo nos altos e media nos baixos
## VII  - media nos altos maximo nos baixos

I <- c(max(matriz[, altos]), max(matriz[, baixo]))
II <- c(max(matriz[, altos]), min(matriz[, baixo]))
III <- c(min(matriz[, altos]), max(matriz[, baixo]))
IV <- c(min(matriz[, altos]), min(matriz[, baixo]))
V <- c(mean(matriz[, altos]), mean(matriz[, baixo]))
VI <- c(max(matriz[, altos]), mean(matriz[, baixo]))
VII <- c(mean(matriz[, altos]), max(matriz[, baixo]))

centroides <- rbind(I, II, III, IV, V, VI, VII) # centróides

Resta agora realizar a análise em si, veja que o laço for() repete a análise para vários valores da constante m e de cada uma, armazena as “probabilidades” de associação a cada centróide e o erro do ajustamento!


ms <- seq(1.2, 10, by = .1)
lista <- array(NA, dim = c(nlevels(dados$fatorB), 
                           nrow(centroides), 
                           length(ms)))
erros <- matrix(NA, nrow = length(ms), ncol = 1)
dimnames(lista) <- list(levels(dados$TRAT), 
                        c('I', 'II', 'III', 'IV', 'V', 'VI', 'VII'), 
                        ms)

for(i in 1:length(ms))
  {
    j <- ms[i]
    ## função de cluster 'fuzzy c-means'
    grupamentos <- cmeans(medias, # matriz dos dados
                          centroides, # centróides dos clusters
                          verbose = FALSE, # não imprime histórico
                          dist = 'euclidean', # distância empregada
                          m = j, # constante de desfuzzyficação
                          method = 'cmeans') # método
    lista[,,i] <- grupamentos$membership
    erros[i,] <- grupamentos$withinerror
  }

Pronto… feito, quer dizer, quase! Realizamos portanto a análise para os vários valores de m resta agora “escolher” um deles e partir para o resto! Até a próxima…

PCA de três dimensões

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

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

Um pouco de teoria:

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

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

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

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

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

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

require(PTAk) # pacote analises pscicrometricas

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

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

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

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

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

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

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

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

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

UFA!

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

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

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

ajuste.PCA <- otimim.PCAn(Z)

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

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

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

Voi lá…, Até mais pessoal!

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

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

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.

Funções para análise de experimentos balanceados

ExpDes é um pacote não oficial do R desenvolvido por Denismar Alves Nogueira, Eric Batista Ferreira e Portya Piscitelli Cavalcanti. O pacote contém funções para análise de dados de experimentos balanceados. As funções são de fácil aplicação e apresentam como resultado o quadro de análise de variância, teste de normalidade para os resíduos, o teste de médias para comparar níveis dos fatores (quando qualitativos) e ajuste de modelos de regressão polinomial (quando quantitativo).

As funções disponíveis no pacote ExpDes (nomes em português) podem ser aplicadas aos seguintes tipos de delineamentos experimentais, sempre considerando fatores de efeito fixo:

função descrição
dic() Delineamento Inteiramente Casualizado balanceado com um só fator.
dbc() Delineamento em Blocos Casualizados balanceado com um só fator.
dql() Delineamento em Quadrado Latino balanceado com um só fator.
fat2.dic() Delineamento Inteiramente Casualizado balanceado em fatorial duplo.
fat2.dbc() Delineamento em Blocos Casualizados balanceado em fatorial duplo.
fat2.ad.dic() Delineamento Inteiramente Casualizado balanceado em fatorial duplo com um tratamento adicional.
fat2.ad.dbc() Delineamento em Blocos Casualizados balanceado em fatorial duplo com um tratamento adicional.
fat3.dic() Delineamento Inteiramente Casualizado balanceado em fatorial triplo.
fat3.dbc() Delineamento em Blocos Casualizados balanceado em fatorial triplo.
fat3.ad.dic() Delineamento Inteiramente Casualizado balanceado em fatorial triplo com um tratamento adicional.
fat3.ad.dbc() Delineamento em Blocos Casualizados balanceado em fatorial triplo com um tratamento adicional.
psub2.dic() Delineamento Inteiramente Casualizado balanceado em esquema de parcelas subdivididas.
psub2.dbc() Delineamento em Blocos Casualizados balanceado em esquema de parcelas subdivididas.

Para casos experimentais onde houve perda de observações (desbalanceados), experimentos com blocos incompletos, experimentos em parcela subsubvividida, experimentos com mais de um tratamento adicional, experimentos com emprego de regressão não linear, experimentos com fatores de efeito aleatório, experimentos com emprego de confundimento entre fatores (fatoriais fracionários), experimentos com resposta não normal (e.g. proporções, contagens) pode-se usar as funções lm(), aov(), nls(), e nlme::lme(), glm() do R, dentre outras. No entanto, o usuário deve ter um conhecimento mais aprofundado sobre o tipo de experimento em questão para evitar conclusões equivocadas dos seus resultados. Esses tipos de experimentos exigem mais conhecimento de modelos lineares, álgebra de matrizes, estimação e inferência.

O pacote está disponível nas versões português e inglês. Confesso que só consegui instalar a versão em inglês mas sei que pessoas que instalaram a versão português com sucesso. Sempre ao usar um pacote, faça a citação do mesmo. Use a função citation("ExpDes"). No link está disponível o manual em português do pacote.