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.
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".
rm(list = ls(all = TRUE)); ls()
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