Senhores,<div><br></div><div>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).</div><div><br></div><div>
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).</div>
<div><br></div><div>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.</div>
<div><br></div><div>Meu ponto de partida é o post do Walmes no wiki do LEG (<a href="http://www.leg.ufpr.br/doku.php/ridiculas#analise_intrabloco_e_interbloco_de_experimento_em_blocos_incompletos_tipo_i">http://www.leg.ufpr.br/doku.php/ridiculas#analise_intrabloco_e_interbloco_de_experimento_em_blocos_incompletos_tipo_i</a>). Todavia, minha implementação não "concorda" com o que lá é mostrado.</div>
<div><br></div><div>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.</div>
<div><br></div><div>Assim, gostaria da ajuda de vocês para aprender um modo menos custoso de reproduzir o jeito original!</div><div><br></div><div>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".</div>
<div><br></div><div>att,</div><div>FH</div><div><br></div><div>Segue o CMR: </div><div><br></div><div><div>rm(list = ls(all = TRUE)); ls()</div><div><br></div><div>download <- read.table('<a href="http://dl.dropbox.com/u/38195533/dados_CC.txt">http://dl.dropbox.com/u/38195533/dados_CC.txt</a>',</div>
<div> header = TRUE, # com cabecalho</div><div> sep = '\t', # separador de celulas <TAB></div><div> dec = ',', # separador de decimal <, (virgula)></div>
<div> na.string = '.') # indicador de omissao</div><div><br></div><div>## lendo e transformando variaveis em fator</div><div>dados <- transform(download, </div><div> trat = factor(trat), # transforma em fator</div>
<div> rep = factor(rep), # idem</div><div> bloco = factor(bloco)) # ...</div><div><br></div><div>str(dados) # estrutura da planilha</div><div><br></div><div>## fator bloco dentro de repeticoes</div>
<div>dados$blocoDrep <- factor(with(dados, paste(rep, bloco, sep = '')))</div><div><br></div><div>## niveis de cada fator **apenas auxiliar**</div><div>r <- nlevels(dados$rep)</div><div>k <- nlevels(dados$bloco)</div>
<div><br></div><div>## array do experimento **generalizacao da matrix do delineamento**</div><div>X <- with(dados, tapply(resp, list(rep, blocoDrep, trat), mean, na.rm = TRUE))</div><div><br></div><div>## totais dos fatores</div>
<div>R <- apply(X, 1, sum, na.rm = TRUE) # repeticoes</div><div>B <- apply(X, 2, sum, na.rm = TRUE) # blocos dentro de repeticoes</div><div>T <- apply(X, 3, sum, na.rm = TRUE) # tratamentos</div><div><br></div><div>
## somas de quadrados do 'jetinho' Pimentel Gomes</div><div>Sx <- sum(X, na.rm = TRUE)</div><div>Sx2 <- sum(X^2, na.rm = TRUE)</div><div>C <- Sx^2 / sum(!<a href="http://is.na">is.na</a>(X))</div><div>SQtotal <- Sx2 - C</div>
<div>SQr <- ( 1 / k^2 ) * sum( R^2 ) - C</div><div>SQt <- ( 1 / r ) * sum( T^2 ) - C</div><div>SQb.unadj <- ( 1 / k ) * sum( B^2 ) - C</div><div>SQe.DBC <- SQtotal - (SQr + SQt)</div><div><br></div><div>## totais dos blocos onde os tratamentos ocorrem... ver Cochran e Cox (cap 10, pag 274)</div>
<div>X.tb <- apply(X, c(2, 3), function(x) sum(!<a href="http://is.na">is.na</a>(x)))</div><div>Bt <- t(X.tb) %*% B</div><div>W <- k * T - (k + 1) * Bt + sum(T)</div><div><br></div><div>## outras somas de quadrados... blocos ajustados e erros</div>
<div>SQb.adj <- sum(W^2) / ( k^3 * ( k + 1 ) )</div><div>SQe.intra <- SQe.DBC - SQb.adj</div><div>Eb <- SQb.adj / (r * ( k - 1 ))</div><div>Ee <- SQe.intra / ( ( k - 1 ) * ( r * k - k - 1 ) )</div><div><br></div>
<div>## coeficiente 'mu' determina se existe ou nao ajuste das medias</div><div>if( Eb <= Ee ) mu <- 0 else mu <- ( Eb - Ee ) / ( k * ( r - 1 ) * Eb )</div><div><br></div><div>Eel <- Ee * ( 1 + ( r * k * mu ) / ( k + 1 ))</div>
<div>EDBC <- SQe.DBC / ( ( r - 1 ) * ( k^2 - 1 ) )</div><div><br></div><div>## soma de quadrados de tratamentos ajustados **PROBLEMA**</div><div>SQt.adj <- sum((T + mu * W - mean(T + mu * W))^2) / r</div><div><br></div>
<div>## montando quadro da ANOVA</div><div>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(!<a href="http://is.na">is.na</a>(X)) - 1 ))</div>
<div>SSs <- c(SQr, SQt, SQt.adj, SQb.adj, SQe.intra, SQe.DBC, Eel * (( k - 1 ) * ( r * k - k - 1 )), SQtotal)</div><div>table <- cbind(DofF, SSs, SSs / DofF)</div><div>dimnames(table) <- list(c('Rep', 'Trat. unadj.', ' Adj.', 'Bloco Adj.', 'Error Intra-bloco', ' DBC', ' Effective', 'Total'), c('Df', 'Sum Sq', 'Mean Sq'))</div>
<div><br></div><div>## teste F para tratamentos ajustados</div><div>F.trat <- table[3, 3] / table[7, 3]</div><div>p.trat <- pf(F.trat, table[3, 1], table[7, 1], lower.tail = FALSE)</div><div><br></div><div>## ajuste via aov</div>
<div>ajuste <- aov(terms(resp ~ rep/bloco + trat, keep = TRUE),</div><div> contrasts = list(rep = contr.sum, bloco = contr.sum, trat = contr.sum),</div><div> data = dados)</div><div><br></div>
<div>## comparando resultados</div><div><br></div><div>## ANOVAS</div><div>table # anova 'made by hands'</div><div>c('F value' = F.trat, 'Pr(>F)' = p.trat) # teste F 'made by hands'</div>
<div>summary(ajuste) # anova via aov()</div><div><br></div><div>## sugestao/implementacao do Walmes (Wiki do LEG)</div><div>car::Anova(ajuste, type = 'III') </div><div>drop1(ajuste, scope = .~., test = 'F')</div>
<div><br></div><div>## medias de tratamentos ajustadas</div><div>data.frame(Unadj. = T / r,</div><div> Adj. = (T + mu * W) / r,</div><div> ## sugestao/implementacao do Walmes (Wiki do LEG)</div><div> Outro = coef(ajuste)[1] + ajuste$contrasts$trat %*% coef(ajuste)[ajuste$assign == 3])</div>
<div><br></div><div>## \EOF</div></div>