[R-br] Fatorial tiplo com dados adicionais

[walmes .]

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.
​
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20150510/2c9e81b5/attachment.html>