[R-br] Fatorial tiplo com dados adicionais

Alisson Lucrécio alisson.lucrecio em ifgoiano.edu.br
Segunda Maio 11 00:13:04 BRT 2015


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


Mais detalhes sobre a lista de discussão R-br