[R-br] Fatorial tiplo com dados adicionais
Alisson Lucrécio
alisson.lucrecio em ifgoiano.edu.br
Quarta Maio 27 21:57:40 BRT 2015
Caro Walmes,
Boa noite.
Quando tenho dados adcionais é possível fazer teste de comparação de
médias, e estudar a interação? Segue adaptação da rotina que passou para os
meus dados.
Obrigado.
Alisson Lucrécio da Costa
#########################################
full_data <- structure(list(Solo = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("LHd", "LVd1", "LVd2",
"OX", "RQ", "SXd"), class = "factor"), Sol = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Água",
"PO4", "SO4"), class = "factor"), Conc = structure(c(1L, 1L,
1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 2L, 2L, 2L, 3L, 3L, 3L,
4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 2L,
2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L,
3L, 4L, 4L, 4L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L,
2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 2L, 2L, 2L, 3L, 3L, 3L, 4L,
4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 2L, 2L,
2L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L,
4L, 4L, 4L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("0",
"0.1", "0.25", "0.5"), class = "factor"), Carbono = c(27.087,
30.352, 29.322, 242, 236.489, 222.298, 432, 438.257, 447.383,
627.138, 631.305, 623.106, 136.73, 75.979, 70.907, 102.099, 109.085,
103.674, 123.428, 128.664, 128.844, 39.875, 27.29, 27.084, 75.825,
79.792, 80.901, 103.56, 107.515, 102.137, 138.102, 136, 141.998,
40.357, 27.128, 30.341, 39.267, 42, 40.453, 46.644, 46.076, 45.548,
74.694, 76.737, 75.685, 395, 403.542, 395.252, 527, 539.921,
530.46, 622.059, 631.452, 629.786, 259.416, 201.837, 200.237,
269.717, 270.958, 269.994, 292, 286.941, 290.546, 72.797, 48.767,
47.596, 187, 181.801, 192.776, 234.527, 232, 229.607, 240, 229.223,
236.598, 79.515, 66.552, 61.401, 92.569, 96.655, 93.394, 104.37,
107.796, 106.256, 75.803, 70.931, 68.759, 391, 414.772, 431.379,
602.842, 615.489, 600, 787.489, 796.342, 798, 110, 134.839, 88.296,
104.46, 100.025, 102.778, 133.712, 136.63, 137.27, 153.782, 155.091,
152.531, 743, 763.369, 752.123, 986.662, 991.683, 1003, 1343.37,
1293, 1330.571, 280, 274.604, 279.925, 314.887, 319.158, 325.88,
350.673, 357.058, 346.08)), .Names = c("Solo", "Sol", "Conc",
"Carbono"), row.names = c(NA, -126L), class = "data.frame")
str(full_data)
xyplot(Carbono ~ Conc|Solo, data=full_data, group = Sol, type=c("p","a"))
ftable(xtabs(~Solo+Sol+Conc, data=full_data))
lm0 <- lm(Carbono ~ Solo*Sol*Conc, data = full_data)
coef(lm0)
anova(lm0)
estm <- names(coef(lm0))[!is.na(coef(lm0))]
X <- model.matrix(lm0)[, estm]
lm1 <- lm(Carbono~0+X, data=full_data)
names(lm1$coefficients) <- colnames(X)
summary(lm1)
c(deviance(lm0), deviance(lm1))
M1 <- by(data=X, INDICES=full_data$Solo, FUN=as.matrix)
M2 <- by(data=X, INDICES=full_data$Sol, FUN=as.matrix)
M3 <- by(data=X, INDICES=full_data$Conc, FUN=as.matrix)
K1 <- t(sapply(M1, FUN=colMeans))
K2 <- t(sapply(M2, FUN=colMeans))
K3 <- t(sapply(M3, FUN=colMeans))
str(K1)
str(K2)
str(K2)
MASS::fractions(t(K1))
MASS::fractions(t(K2))
MASS::fractions(t(K3))
K1%*%coef(lm1)
K2%*%coef(lm1)
K3%*%coef(lm1)
ddply(full_data, "Solo", summarise, Carbono = mean(Carbono))
ddply(full_data, "Sol", summarise, Carbono = mean(Carbono))
ddply(full_data, "Conc", summarise, Carbono = mean(Carbono))
G1 <- apc(K1)
G2 <- apc(K2)
G3 <- apc(K3)
summary(glht(lm1, linfct=G1), test=adjusted(type="fdr"))
summary(glht(lm1, linfct=G2), test=adjusted(type="fdr"))
summary(glht(lm1, linfct=G3), test=adjusted(type="fdr"))
################
2015-05-11 0:13 GMT-03:00 Alisson Lucrécio <alisson.lucrecio em ifgoiano.edu.br
>:
> Obrigado, Walmes.
>
> Vou analisar a rotina.
>
> Att.
>
> 2015-05-10 23:10 GMT-03:00 Walmes Zeviani 2 [via R-br] <
> ml-node+s2285057n4664466h82 em n4.nabble.com>:
>
>> 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.
>>
>>
>> _______________________________________________
>> R-br mailing list
>> [hidden email] <http:///user/SendEmail.jtp?type=node&node=4664466&i=0>
>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>> código mínimo reproduzível.
>>
>> ------------------------------
>> If you reply to this email, your message will be added to the
>> discussion below:
>>
>> http://r-br.2285057.n4.nabble.com/R-br-Fatorial-tiplo-com-dados-adicionais-tp4664459p4664466.html
>> To unsubscribe from R-br, click here
>> <http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3357982&code=YWxpc3Nvbi5sdWNyZWNpb0BpZmdvaWFuby5lZHUuYnJ8MzM1Nzk4MnwtNTI1NDI5NDE3>
>> .
>> NAML
>> <http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>>
>
>
>
> --
> Alisson Lucrecio da Costa
>
--
Alisson Lucrecio da Costa
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20150527/c4335ce1/attachment.html>
Mais detalhes sobre a lista de discussão R-br