<div dir="ltr">Caro Walmes,<div><br></div><div>Boa noite.</div><div><br></div><div>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.</div><div><br></div><div>Obrigado.</div><div><br></div><div>Alisson Lucrécio da Costa</div><div><br></div><div>#########################################</div><div><br></div><div>full_data <- structure(list(Solo = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 4L,</div><span class="" style="border-collapse:separate;color:rgb(0,0,0);font-family:'Ubuntu Mono';font-size:14px;line-height:16px;white-space:pre-wrap;border-spacing:0px;background-color:rgb(225,226,229)"><pre tabindex="0" class="" id="rstudio_console_output" style="font-family:'Ubuntu Mono';outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px;line-height:1.2;font-size:10.4pt!important;white-space:pre-wrap!important">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")</pre></span><div><br></div><div>str(full_data)<br></div><div><br></div><div><div>xyplot(Carbono ~ Conc|Solo, data=full_data, group = Sol, type=c("p","a"))</div><div><br></div><div>ftable(xtabs(~Solo+Sol+Conc, data=full_data))</div><div><br></div><div>lm0 <- lm(Carbono ~ Solo*Sol*Conc, data = full_data)</div><div><br></div><div>coef(lm0)</div><div><br></div><div>anova(lm0)</div><div><br></div><div>estm <- names(coef(lm0))[!<a href="http://is.na">is.na</a>(coef(lm0))]</div><div><br></div><div>X <- model.matrix(lm0)[, estm]</div><div><br></div><div>lm1 <- lm(Carbono~0+X, data=full_data)</div><div><br></div><div>names(lm1$coefficients) <- colnames(X)</div><div><br></div><div>summary(lm1)</div><div><br></div><div>c(deviance(lm0), deviance(lm1))</div><div><br></div><div>M1 <- by(data=X, INDICES=full_data$Solo, FUN=as.matrix)</div><div>M2 <- by(data=X, INDICES=full_data$Sol, FUN=as.matrix)</div><div>M3 <- by(data=X, INDICES=full_data$Conc, FUN=as.matrix)</div><div><br></div><div>K1 <- t(sapply(M1, FUN=colMeans))</div><div>K2 <- t(sapply(M2, FUN=colMeans))</div><div>K3 <- t(sapply(M3, FUN=colMeans))</div><div><br></div><div>str(K1)</div><div>str(K2)</div><div>str(K2)</div><div><br></div><div>MASS::fractions(t(K1))</div><div>MASS::fractions(t(K2))</div><div>MASS::fractions(t(K3))</div><div><br></div><div>K1%*%coef(lm1)</div><div>K2%*%coef(lm1)</div><div>K3%*%coef(lm1)</div><div><br></div><div>ddply(full_data, "Solo", summarise, Carbono = mean(Carbono))</div><div>ddply(full_data, "Sol", summarise, Carbono = mean(Carbono))</div><div>ddply(full_data, "Conc", summarise, Carbono = mean(Carbono))</div><div><br></div><div>G1 <- apc(K1)</div><div>G2 <- apc(K2)</div><div>G3 <- apc(K3)</div><div><br></div><div>summary(glht(lm1, linfct=G1), test=adjusted(type="fdr"))</div><div>summary(glht(lm1, linfct=G2), test=adjusted(type="fdr"))</div><div>summary(glht(lm1, linfct=G3), test=adjusted(type="fdr"))</div></div><div><br></div><div>################</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-05-11 0:13 GMT-03:00 Alisson Lucrécio <span dir="ltr"><<a href="mailto:alisson.lucrecio@ifgoiano.edu.br" target="_blank">alisson.lucrecio@ifgoiano.edu.br</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Obrigado, Walmes.<div><br></div><div>Vou analisar a rotina.</div><div><br></div><div>Att.</div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-05-10 23:10 GMT-03:00 Walmes Zeviani 2 [via R-br] <span dir="ltr"><<a href="mailto:ml-node+s2285057n4664466h82@n4.nabble.com" target="_blank">ml-node+s2285057n4664466h82@n4.nabble.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><div><div>
<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" rel="nofollow" link="external" target="_blank">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" rel="nofollow" link="external" target="_blank">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" rel="nofollow" link="external" target="_blank">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" rel="nofollow" link="external" target="_blank">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>
<br></div></div></div></div><span class="">_______________________________________________
<br>R-br mailing list
<br><a href="http:///user/SendEmail.jtp?type=node&node=4664466&i=0" rel="nofollow" link="external" target="_blank">[hidden email]</a>
<br><a href="https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br" rel="nofollow" link="external" target="_blank">https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br</a><br>Leia o guia de postagem (<a href="http://www.leg.ufpr.br/r-br-guia" rel="nofollow" link="external" target="_blank">http://www.leg.ufpr.br/r-br-guia</a>) e forneça código mínimo reproduzível.
<br>
<br>
<hr noshade size="1" color="#cccccc">
<div style="color:#444;font:12px tahoma,geneva,helvetica,arial,sans-serif">
<div style="font-weight:bold">If you reply to this email, your message will be added to the discussion below:</div>
<a href="http://r-br.2285057.n4.nabble.com/R-br-Fatorial-tiplo-com-dados-adicionais-tp4664459p4664466.html" target="_blank">http://r-br.2285057.n4.nabble.com/R-br-Fatorial-tiplo-com-dados-adicionais-tp4664459p4664466.html</a>
</div>
<div style="color:#666;font:11px tahoma,geneva,helvetica,arial,sans-serif;margin-top:.4em;line-height:1.5em">
To unsubscribe from R-br, <a href="http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3357982&code=YWxpc3Nvbi5sdWNyZWNpb0BpZmdvaWFuby5lZHUuYnJ8MzM1Nzk4MnwtNTI1NDI5NDE3" target="_blank">click here</a>.<br>
<a href="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" rel="nofollow" style="font:9px serif" target="_blank">NAML</a>
</div></span></blockquote></div><span class="HOEnZb"><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div><div dir="ltr"><span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:13px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;display:inline!important;float:none;background-color:rgb(255,255,255)">Alisson Lucrecio da Costa</span></div></div>
</font></span></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature"><div dir="ltr"><span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:13px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;display:inline!important;float:none;background-color:rgb(255,255,255)">Alisson Lucrecio da Costa</span></div></div>
</div>