[R-br] Lattice Design

[FHRB Toledo]

Senhores,

Estou com uma dúvida com relação a reprodução de uma análise de um
experimento em delineamento em lattice quadrado balanceado, tipo I, segundo
Pimentel Gomes (cap 10).

Escreve esse email baseado tanto no livro do Pimentel Gomes, como também no
Cochran e Cox, cap 10, pag 274. Para quem interessar segue também um
exemplo (teste) do livro Ramalho et al., Experimentação em Genética ... cap
11, pág 171 (mas não é estelionato, rsrs).

Pois bem, como disse, trata-se da análise de blocos incompletos em lattice
quadrado balanceado. A pergunta é: como obter a soma de quadrados de
tratamentos ajustadas aos efeitos de blocos? E ainda, de que forma obter
também as médias de tratamentos livre dos efeitos de blocos.

Meu ponto de partida é o post do Walmes no wiki do LEG (
http://www.leg.ufpr.br/doku.php/ridiculas#analise_intrabloco_e_interbloco_de_experimento_em_blocos_incompletos_tipo_i).
Todavia, minha implementação não "concorda" com o que lá é mostrado.

Peço desculpas pela forma que procedi, mas a título didático me prendi
exatamente ao modo como é originalmente apresentado em Cochran e Cox (pag
274) e seguindo o que lá é feito consegui reproduzir ipsis literis a
análise.

Assim, gostaria da ajuda de vocês para aprender um modo menos custoso de
reproduzir o jeito original!

Desde já agradeço qualquer ajuda/sugestão com relação a dúvida! Caso alguém
se interesse disponibilizo outro conjunto de dados, da mesma natureza no
mesmo lik, com o nome de "dados_magno.txt".

att,
FH

Segue o CMR:

rm(list = ls(all = TRUE)); ls()

download <- read.table('http://dl.dropbox.com/u/38195533/dados_CC.txt',
                       header = TRUE, # com cabecalho
                       sep = '\t', # separador de celulas <TAB>
                       dec = ',', # separador de decimal <, (virgula)>
                       na.string = '.') # indicador de omissao

## lendo e transformando variaveis em fator
dados <- transform(download,
                   trat = factor(trat), # transforma em fator
                   rep = factor(rep), # idem
                   bloco = factor(bloco)) # ...

str(dados) # estrutura da planilha

## fator bloco dentro de repeticoes
dados$blocoDrep <- factor(with(dados, paste(rep, bloco, sep = '')))

## niveis de cada fator **apenas auxiliar**
r <- nlevels(dados$rep)
k <- nlevels(dados$bloco)

## array do experimento **generalizacao da matrix do delineamento**
X <- with(dados, tapply(resp, list(rep, blocoDrep, trat), mean, na.rm =
TRUE))

## totais dos fatores
R <- apply(X, 1, sum, na.rm = TRUE) # repeticoes
B <- apply(X, 2, sum, na.rm = TRUE) # blocos dentro de repeticoes
T <- apply(X, 3, sum, na.rm = TRUE) # tratamentos

## somas de quadrados do 'jetinho' Pimentel Gomes
Sx <- sum(X, na.rm = TRUE)
Sx2 <- sum(X^2, na.rm = TRUE)
C <- Sx^2 / sum(!is.na(X))
SQtotal <- Sx2 - C
SQr <- ( 1 / k^2 ) * sum( R^2 ) - C
SQt <- ( 1 / r ) * sum( T^2 ) - C
SQb.unadj <- ( 1 / k ) * sum( B^2 ) - C
SQe.DBC <- SQtotal - (SQr + SQt)

## totais dos blocos onde os tratamentos ocorrem... ver Cochran e Cox (cap
10, pag 274)
X.tb <- apply(X, c(2, 3), function(x) sum(!is.na(x)))
Bt <- t(X.tb) %*% B
W <- k * T - (k + 1) * Bt + sum(T)

## outras somas de quadrados... blocos ajustados e erros
SQb.adj <- sum(W^2) / ( k^3 * ( k + 1 ) )
SQe.intra <- SQe.DBC - SQb.adj
Eb <- SQb.adj / (r * ( k - 1 ))
Ee <- SQe.intra / ( ( k - 1 ) * ( r * k - k - 1 ) )

## coeficiente 'mu' determina se existe ou nao ajuste das medias
if( Eb <= Ee ) mu <- 0 else mu <- ( Eb - Ee ) / ( k * ( r - 1 ) * Eb )

Eel <- Ee * ( 1 + ( r * k * mu ) / ( k + 1 ))
EDBC <- SQe.DBC / ( ( r - 1 ) * ( k^2 - 1 ) )

## soma de quadrados de tratamentos ajustados **PROBLEMA**
SQt.adj <- sum((T + mu * W - mean(T + mu * W))^2) / r

## montando quadro da ANOVA
DofF <- c( (r - 1 ), ( k^2 - 1 ), ( k^2 - 1 ), ( r * ( k - 1 ) ), ( ( k - 1
) * ( r * k - k - 1 ) ), ( ( r - 1 ) * ( k^2 - 1 ) ), ( ( k - 1 ) * ( r * k
- k - 1 ) ), ( sum(!is.na(X)) - 1 ))
SSs <- c(SQr, SQt, SQt.adj, SQb.adj, SQe.intra, SQe.DBC, Eel * (( k - 1 ) *
( r * k - k - 1 )), SQtotal)
table <- cbind(DofF, SSs, SSs / DofF)
dimnames(table) <- list(c('Rep', 'Trat. unadj.', '      Adj.', 'Bloco
Adj.', 'Error Intra-bloco', '      DBC', '      Effective', 'Total'),
c('Df', 'Sum Sq', 'Mean Sq'))

## teste F para tratamentos ajustados
F.trat <- table[3, 3] / table[7, 3]
p.trat <- pf(F.trat, table[3, 1], table[7, 1], lower.tail = FALSE)

## ajuste via aov
ajuste <- aov(terms(resp ~ rep/bloco + trat, keep = TRUE),
              contrasts = list(rep = contr.sum, bloco = contr.sum, trat =
contr.sum),
              data = dados)

## comparando resultados

## ANOVAS
table # anova 'made by hands'
c('F value' = F.trat, 'Pr(>F)' = p.trat) # teste F 'made by hands'
summary(ajuste) # anova via aov()

## sugestao/implementacao do Walmes (Wiki do LEG)
car::Anova(ajuste, type = 'III')
drop1(ajuste, scope = .~., test = 'F')

## medias de tratamentos ajustadas
data.frame(Unadj. = T / r,
           Adj. = (T + mu * W) / r,
           ## sugestao/implementacao do Walmes (Wiki do LEG)
           Outro = coef(ajuste)[1] + ajuste$contrasts$trat %*%
coef(ajuste)[ajuste$assign == 3])

## \EOF
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121115/6e7ae921/attachment.html>