<div dir="ltr"><div><div><div><div><div><div><div><div>Prezados membros da lista,<br><br></div>Alguns dias postei um código não reproduzível sobre modelos para heterocedásticidade.<br></div>Peço desculpas e abaixo vai o código que imagino ser possível reproduzir.<br></div>Os dois principais problemas eram a instalação do pacote mcglm e o conjunto de dados.<br></div>O pacote pode ser instalado facilmente pelo github repository neste endereço<br><br><a href="https://github.com/wbonat/mcglm">https://github.com/wbonat/mcglm</a><br>install_github("wbonat/mcglm", ref = "devel") # Saliento para instalar a versão devel !!<br><br></div>O conjunto de dados está anexado neste e-mail.<br></div>Além disso, o Luiz Leal estava com problemas pra ter acesso e pediu o .tar.gz que também vai em anexo.<br></div>Aproveito para salientar que o mcglm está em desenvolvimento, assim quaisquer dúvidas, criticas e/ou <br>sugestões serão muito bem vindas.<br></div><div><div><br></div><div>All the best!<br></div><div><br># Heteroscedastic regression model -------------------------------------<br># Author: Wagner Hugo Bonat LEG/UFPR -----------------------------------<br># Date: 13/11/2016 -----------------------------------------------------<br><br># Loading extra packages<br>install.packages("devtools")<br>require(devtools)<br><br># Install mcglm package from github repository -------------------------<br>install_github("wbonat/mcglm", ref = "devel")<br>require(mcglm)<br><br><br># Loading data set<br>Fenois <- c(337.311, 344.874, 342.353, 325.546, 333.950, 330.588, 328.067,<br> 328.067, 318.824, 331.429, 333.950, 334.790, 336.471, 338.151,<br> 342.353, 259.160, 252.437, 268.403, 265.882, 266.723, 287.731,<br> 88.571, 88.571, 90.252, 41.513, 52.437, 49.076, 88.571,<br> 88.571, 90.252, 64.202, 60.000, 61.681)<br>Cor <- factor(c(rep("ambar",6), rep("ambar_claro",3), rep("ambar",6),<br> rep("ambar_claro",6),rep("branco",6),<br> rep("extra_ambar_claro",3),rep("branco",3)))<br>dados <- data.frame("Fenois" = Fenois, "Cor" = Cor)<br>boxplot(dados$Fenois ~ dados$Cor)<br>tapply(dados$Fenois, dados$Cor, sd)<br><br># Linear regression model- ---------------------------------------------<br>fit1 <- lm(Fenois ~ Cor, data = dados)<br>anova(fit1)<br>plot(residuals(fit1) ~ fitted(fit1))<br>plot(residuals(fit1) ~ Cor)<br><br># Double Linear regression model ---------------------------------------<br>dados$id <- 1<br>fit2 <- mcglm(c(Fenois ~ Cor), list(mc_dglm(~ Cor, id = "id", data = dados)), <br> covariance = "expm", data = dados)<br>summary(fit2)<br><br>summary(fit1)<br>summary(fit2)<br><br>cbind(coef(fit1), coef(fit2, type = "beta")$Estimates)<br>cbind(sqrt(diag(vcov(fit1))), coef(fit2, type = "beta", std.error = TRUE)$Std.error)<br><br># Example 2 ------------------------------------------------------------<br>dados2 <- read.table("HETE.csv", header = TRUE, sep = ";", dec = ",")<br>with(dados2, boxplot(y ~ x))<br><br># Note that, the A33 has no variance, so we need to remove this level.<br>dados2 <- dados2[which(dados2$x != "A33"),]<br>dados2$x <- droplevels(dados2$x)<br>tapply(dados2$y, dados2$x, sd)<br><br># Linear regression model ----------------------------------------------<br>fit_lm <- lm(y ~ x, data = dados2)<br>summary(fit_lm)<br>plot(residuals(fit_lm) ~ fitted(fit_lm))<br>plot(residuals(fit_lm) ~ dados2$x)<br><br># Double linear regression model ---------------------------------------<br>dados2$id <- 1<br>fit_dlm <- mcglm(c(y ~ x), list(mc_dglm(~ x, id = "id", data = dados2)), <br> covariance = "expm", data = dados2)<br>summary(fit_dlm)<br><br># Comparing estimates and standard errors ------------------------------<br>cbind(coef(fit_lm), coef(fit_dlm, type = "beta")$Estimates)<br>cbind(sqrt(diag(vcov(fit_lm))), coef(fit_dlm, type = "beta", std.error = TRUE)$Std.error)<br><br><div><div><br></div></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">Em 11 de novembro de 2016 17:01, Luiz Leal <span dir="ltr"><<a href="mailto:richfield1974@yahoo.com" target="_blank">richfield1974@yahoo.com</a>></span> escreveu:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div style="color:#000;background-color:#fff;font-family:HelveticaNeue,Helvetica Neue,Helvetica,Arial,Lucida Grande,sans-serif;font-size:16px"><div><span>Prezado Wagner, uso o R no trabalho e devido a restrições técnicas só consigo instalar pacotes a partir do zip. Poderia me mandar o zip deste pacote?</span></div><div><span>Att</span></div><div><span>Luiz</span></div> <div class="m_-4698200378919308510qtdSeparateBR"><br><br></div><div class="m_-4698200378919308510yahoo_quoted" style="display:block"> <div style="font-family:HelveticaNeue,Helvetica Neue,Helvetica,Arial,Lucida Grande,sans-serif;font-size:16px"> <div style="font-family:HelveticaNeue,Helvetica Neue,Helvetica,Arial,Lucida Grande,sans-serif;font-size:16px"><div><div class="h5"> <div dir="ltr"><font face="Arial" size="2"> On Friday, November 4, 2016 5:05 PM, Wagner Bonat <<a href="mailto:wbonat@gmail.com" target="_blank">wbonat@gmail.com</a>> wrote:<br></font></div> <br><br> </div></div><div class="m_-4698200378919308510y_msg_container"><div id="m_-4698200378919308510yiv1337977286"><div><div><div class="h5"><div dir="ltr"><div><div><div>Sim, vc pode usar qualquer dos testes padrões. <br clear="none">A única diferença é que vc deve usar os beta e erros padrões que vem deste modelo.<br clear="none">Não sei se vc pode incluir sua própria matriz de covariancia nas funções do R.<br clear="none">Talvez, sim. Se não tem que implementar.<br clear="none"></div>O Walmes é muito bom nestes testes de comparações múltiplas talvez ele possa ajudar.<br clear="none"></div><div><br clear="none">O que houve q não conseguiu instalar o pacote?<br clear="none"></div><div>Tentou pelo github ou do CRAN?<br clear="none"></div><div><br clear="none"><a rel="nofollow" shape="rect" href="https://cran.r-project.org/web/packages/mcglm/index.html" target="_blank">https://cran.r-project.org/<wbr>web/packages/mcglm/index.html</a><br clear="none"><a rel="nofollow" shape="rect" href="https://github.com/wbonat/mcglm" target="_blank">https://github.com/wbonat/<wbr>mcglm</a><br clear="none"><br clear="none"></div></div><br clear="none"><br clear="none"></div></div></div><div class="m_-4698200378919308510yiv1337977286gmail_extra"><br clear="none"><div class="m_-4698200378919308510yiv1337977286yqt1293506770" id="m_-4698200378919308510yiv1337977286yqt16758"><div class="m_-4698200378919308510yiv1337977286gmail_quote"><div><div class="h5">Em 4 de novembro de 2016 19:46, Luiz Leal <span dir="ltr"><<a rel="nofollow" shape="rect" href="mailto:richfield1974@yahoo.com" target="_blank">richfield1974@yahoo.com</a>></span> escreveu:<br clear="none"></div></div><blockquote class="m_-4698200378919308510yiv1337977286gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div style="color:#000;background-color:#fff;font-family:HelveticaNeue,Helvetica Neue,Helvetica,Arial,Lucida Grande,sans-serif;font-size:16px"><div><div class="h5"><div id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yui_3_16_0_ym19_1_1478284318272_4940"><span id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yui_3_16_0_ym19_1_1478284318272_5020">Wagner, meu interesse é, uma vez identificado que existe diferença entre os tratamentos (considerando que um deles é o controle) utilizar o teste de Dunnett para verificar quais tratamentos diferem do tratamento controle. Como o pressuposto de homogeneidade das variâncias é violado busquei alternativas para "homogeneizar" as variâncias. Posso aplicar esse teste a partir do modelo acima descrito?</span></div><div id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yui_3_16_0_ym19_1_1478284318272_4940"><span>Desde já agradeço</span></div><div id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yui_3_16_0_ym19_1_1478284318272_4940"><span>Luiz</span></div><div id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yui_3_16_0_ym19_1_1478284318272_4940"><span id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yui_3_16_0_ym19_1_1478284318272_5175">PS. Não consegui instalar o pacote</span></div> <div class="m_-4698200378919308510yiv1337977286m_-7769256875373248989qtdSeparateBR"><br clear="none"><br clear="none"></div></div></div><div class="m_-4698200378919308510yiv1337977286m_-7769256875373248989yahoo_quoted" style="display:block"> <div style="font-family:HelveticaNeue,Helvetica Neue,Helvetica,Arial,Lucida Grande,sans-serif;font-size:16px"> <div style="font-family:HelveticaNeue,Helvetica Neue,Helvetica,Arial,Lucida Grande,sans-serif;font-size:16px"><div><div class="h5"><div><div class="m_-4698200378919308510yiv1337977286h5"> <div dir="ltr"><font face="Arial" size="2"> On Friday, November 4, 2016 4:13 PM, Wagner Bonat via R-br <<a rel="nofollow" shape="rect" href="mailto:r-br@listas.c3sl.ufpr.br" target="_blank">r-br@listas.c3sl.ufpr.br</a>> wrote:<br clear="none"></font></div> <br clear="none"><br clear="none"> </div></div></div></div><div class="m_-4698200378919308510yiv1337977286m_-7769256875373248989y_msg_container"><div><div class="m_-4698200378919308510yiv1337977286h5"><div id="m_-4698200378919308510yiv1337977286m_-7769256875373248989yiv3974050584"><div dir="ltr"><div><div class="h5"><div>Caros, <br clear="none"></div><div><br clear="none">Alguém postou esse conjunto de dados com problema de pressupostos, principalmente heterocedasticidade. Agora a pouco veio outro e-mail com um problema similar. Fiz um exemplo um pouco mais detalhado de como isso pode<br clear="none"></div>ser facilmente resolvido e mostrando o efeito disso no modelo.<br clear="none"><br clear="none"># Example 2 ------------------------------ ------------------------------<br clear="none">Fenois = c(337.311, 344.874, 342.353, 325.546, 333.950, 330.588, 328.067, <br clear="none"> 328.067, 318.824, 331.429, 333.950, 334.790, 336.471, 338.151, <br clear="none"> 342.353, 259.160, 252.437, 268.403, 265.882, 266.723, 287.731, <br clear="none"> 88.571, 88.571, 90.252, 41.513, 52.437, 49.076, 88.571, <br clear="none"> 88.571, 90.252, 64.202, 60.000, 61.681)<br clear="none">Cor = factor(c(rep("ambar",6), rep("ambar_claro",3), rep("ambar",6), <br clear="none"> rep("ambar_claro",6),rep(" branco",6), <br clear="none"></div></div> rep("extra_ambar_claro",3), rep("branco",3)))<div><div class="h5"><br clear="none"><br clear="none"># Exploratory analysis<br clear="none">boxplot(Fenois ~ Cor)<br clear="none">tapply(Fenois, Cor, sd)<br clear="none">dados <- data.frame(Fenois, Cor)<br clear="none">dados$id <- 1<br clear="none"><br clear="none"># Fitting <br clear="none">fit1 <- mcglm(c(Fenois ~ Cor), list(mc_id(dados)), data = dados)<br clear="none">fit2 <- mcglm(c(Fenois ~ Cor), list(mc_dglm(~ Cor, id = "id", data = dados)),<br clear="none"> covariance = "expm", data = dados)<br clear="none"># Goodness-of-fit <br clear="none">gof(fit1)<br clear="none">gof(fit2)<br clear="none"><br clear="none"># Comparing estimates and standard errors<br clear="none">coef(fit1, type = "beta", std.error = TRUE)<br clear="none">coef(fit2, type = "beta", std.error = TRUE)<br clear="none"><br clear="none">O interessante é que a estimativa pontual é exatamente a mesma, porém<br clear="none">olha a enorme diferença nos erros padrões dos betas.<br clear="all"><div><div><br clear="none">-- <br clear="none"><div class="m_-4698200378919308510yiv1337977286m_-7769256875373248989yiv3974050584gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>Wagner Hugo Bonat<br clear="none">------------------------------ ------------------------------ ------------------------------ ----<br clear="none">Department of Mathematics and Computer Science (IMADA)<br clear="none">University of Southern Denmark (SDU) and<br clear="none">Laboratório de Estatística e Geoinformação (LEG)<br clear="none">Universidade Federal do Paraná (UFPR)<br clear="none"><br clear="none"></div></div></div></div></div>
</div></div></div></div></div></div><br clear="none"></div></div><div><div class="h5">______________________________ _________________<br clear="none">R-br mailing list<br clear="none"><a rel="nofollow" shape="rect" href="mailto:R-br@listas.c3sl.ufpr.br" target="_blank">R-br@listas.c3sl.ufpr.br</a><br clear="none"><a rel="nofollow" shape="rect" href="https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br" target="_blank">https://listas.inf.ufpr.br/ cgi-bin/mailman/listinfo/r-br</a><br clear="none"></div></div>Leia o guia de postagem (<a rel="nofollow" shape="rect" href="http://www.leg.ufpr.br/r-br-guia" target="_blank">http://www.leg.ufpr.br/r-br- guia</a>) e forne� c�igo m�imo reproduz�el.<br clear="none"><br clear="none"></div> </div> </div> </div></div></div></blockquote></div></div><span class=""><br clear="none"><br clear="all"><br clear="none">-- <br clear="none"><div class="m_-4698200378919308510yiv1337977286gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>Wagner Hugo Bonat<br clear="none">------------------------------<wbr>------------------------------<wbr>------------------------------<wbr>----<br clear="none">Department of Mathematics and Computer Science (IMADA)<br clear="none">University of Southern Denmark (SDU) and<br clear="none">Laboratório de Estatística e Geoinformação (LEG)<br clear="none">Universidade Federal do Paraná (UFPR)<br clear="none"><br clear="none"></div></div></div></div></div>
</span></div></div></div><br><br></div> </div> </div> </div></div></div></blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>Wagner Hugo Bonat<br>----------------------------------------------------------------------------------------------<br>Department of Mathematics and Computer Science (IMADA)<br>University of Southern Denmark (SDU) and<br>Laboratório de Estatística e Geoinformação (LEG)<br>Universidade Federal do Paraná (UFPR)<br><br></div></div></div></div></div>
</div>