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