<div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Segue uma rotina que pode ser útil ou servir de inspiração.<br><br><span style="font-family:monospace,monospace">##-----------------------------------------------------------------------------<br><br>require(doBy)<br>require(multcomp)<br><br>str(fat_data)<br>str(adi_data)<br><br>adi_data$Conc <- factor(0, levels=c(0, levels(fat_data$Conc)))<br><br>da <- merge(fat_data, adi_data, all=TRUE)<br>str(da)<br><br>da$Conc <- factor(da$Conc,<br>                  levels=sort(as.numeric(levels(da$Conc))))<br>da$Solo <- factor(da$Solo, levels=c("E","A","B","C","D"))<br>da$Sol <- factor(da$Sol, levels=c("Água","PO4","SO4"))<br>str(da)<br><br>ftable(xtabs(~Solo+Sol+Conc, data=da))<br><br>## Especificação do modelo completo.<br>m0 <- lm(Rep~Solo*Sol*Conc, data=da)<br><br>## Efeitos não estimáveis devido ausência de celas.<br>coef(m0)<br><br>## Quadro de anova.<br>anova(m0)<br><br>##--------------------------------------------<br>## TRUQUE: criar um conjunto de dados artificiais mas que seja um<br>## fatorial completo.<br><br>## Dados artificiais.<br>fac <- c("Solo","Sol","Conc")<br>L <- lapply(fac, function(x) levels(da$x))<br><br>m0$xlevels<br><br>## Dados artificiais, possui todas as celas.<br>db <- do.call(what=expand.grid, args=m0$xlevels)<br>db$Rep <- runif(nrow(db))<br><br>## Fatorial completo com um registro por cela.<br>ftable(xtabs(~Solo+Sol+Conc, data=db))<br><br>## Com este, ajustar um modelo de mentira apenas para que se possa obter<br>## a matriz de coeficientes para se chegar as médias marginais.<br>## mb <- lm(formula(m0), data=db)<br>mb <- update(m0, data=db)<br>anova(mb)<br><br>## Ajustou sem restar graus de liberdade pois não colocou-se<br>## repetições. Não é um problema pois quer-se apenas tirar proveito da<br>## estrutura completa.<br><br>sum(<a href="http://is.na">is.na</a>(coef(mb))) ## Todos os efeitos estimados.<br><br>## Matriz para médias ajustadas de Solo.<br>t(LSmatrix(mb, effect="Solo"))<br><br>## Cuidado! Essa matriz não contém os pesos corretos. É muito útil em<br>## experimentos fatorias completos, mas nos incompletos deve ser usada<br>## com o devido cuidado.<br><br>## Nessa matriz, remover as colunas dos efeitos não estimáveis. Guardar<br>## o nome dos efeitos estimados.<br>estm <- names(coef(m0))[!<a href="http://is.na">is.na</a>(coef(m0))]<br><br>## Matriz do modelo com colunas correspondentes à efeitos estimáveis.<br>X <- model.matrix(m0)[, estm]<br><br>## Reajuste do modelo sem uso de fórmula, mas com o uso da matriz do<br>## modelo contendo apenas colunas de efeitos estimáveis.<br>m1 <- lm(Rep~0+X, data=da)<br><br>c(deviance(m0), deviance(m1)) ## São o mesmo modelo.<br><br>## Partindo a matriz do modelo. Havendo desbalanceamento (caselas<br>## presentes com frequência não igual), esse passo precisa ser revisto e<br>## adaptado.<br><br>M <- by(data=X, INDICES=da$Solo, FUN=as.matrix)<br>## M <- by(data=X, INDICES=da$Sol, FUN=as.matrix)<br>## M <- by(data=X, INDICES=da$Conc, FUN=as.matrix)<br><br>K <- t(sapply(M, FUN=colMeans))<br>str(K)<br><br>## "Pesos". Os pesos estão corretos, diferente daquela retornada pela<br>## LSmatrix().<br>MASS::fractions(t(K))<br><br>## Médias ajustadas.<br>K%*%coef(m1)<br><br>## Médias amostrais coincidem com as médias matriciais.<br>with(da, tapply(Rep, Solo, FUN=mean))<br>## with(da, tapply(Rep, Sol, FUN=mean))<br>## with(da, tapply(Rep, Conc, FUN=mean))<br><br>## Médias ajustadas são sempre bem vindas pela maioria das pessoas mas<br>## elas raramente fazem sentido. Por exemplo, fazem sentido se os<br>## efeitos marginalizados são nulos (termo pode ser removido do modelo)<br>## ou talvez considerados como aleatórios (esperança 0).<br><br>##--------------------------------------------<br>## Contrastes entre níveis de Solo.<br><br>require(wzRfun)<br><br>## No caso de não desejar instalar o wzRfun, apenas copie o código da<br>## função apc() para uma sessão R. Acesso pelo link:<br>## <a href="https://raw.githubusercontent.com/walmes/wzRfun/master/R/apc.R">https://raw.githubusercontent.com/walmes/wzRfun/master/R/apc.R</a><br><br>source("<a href="https://raw.githubusercontent.com/walmes/wzRfun/master/R/apc.R">https://raw.githubusercontent.com/walmes/wzRfun/master/R/apc.R</a>")<br><br>G <- apc(K)<br><br>summary(glht(m1, linfct=G), test=adjusted(type="fdr"))<br>## summary(glht(m1, linfct=G), test=adjusted(type="bonferroni"))<br><br>##-----------------------------------------------------------------------------</span><br><br>À disposição.<br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Walmes.<br></div>​</div>