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