[R-br] Heterocedasticidade e teste de média
Wagner Bonat
wbonat em gmail.com
Domingo Novembro 13 15:40:04 BRST 2016
Prezados membros da lista,
Alguns dias postei um código não reproduzível sobre modelos para
heterocedásticidade.
Peço desculpas e abaixo vai o código que imagino ser possível reproduzir.
Os dois principais problemas eram a instalação do pacote mcglm e o conjunto
de dados.
O pacote pode ser instalado facilmente pelo github repository neste endereço
https://github.com/wbonat/mcglm
install_github("wbonat/mcglm", ref = "devel") # Saliento para instalar a
versão devel !!
O conjunto de dados está anexado neste e-mail.
Além disso, o Luiz Leal estava com problemas pra ter acesso e pediu o
.tar.gz que também vai em anexo.
Aproveito para salientar que o mcglm está em desenvolvimento, assim
quaisquer dúvidas, criticas e/ou
sugestões serão muito bem vindas.
All the best!
# Heteroscedastic regression model -------------------------------------
# Author: Wagner Hugo Bonat LEG/UFPR -----------------------------------
# Date: 13/11/2016 -----------------------------------------------------
# Loading extra packages
install.packages("devtools")
require(devtools)
# Install mcglm package from github repository -------------------------
install_github("wbonat/mcglm", ref = "devel")
require(mcglm)
# Loading data set
Fenois <- c(337.311, 344.874, 342.353, 325.546, 333.950, 330.588, 328.067,
328.067, 318.824, 331.429, 333.950, 334.790, 336.471, 338.151,
342.353, 259.160, 252.437, 268.403, 265.882, 266.723, 287.731,
88.571, 88.571, 90.252, 41.513, 52.437, 49.076, 88.571,
88.571, 90.252, 64.202, 60.000, 61.681)
Cor <- factor(c(rep("ambar",6), rep("ambar_claro",3), rep("ambar",6),
rep("ambar_claro",6),rep("branco",6),
rep("extra_ambar_claro",3),rep("branco",3)))
dados <- data.frame("Fenois" = Fenois, "Cor" = Cor)
boxplot(dados$Fenois ~ dados$Cor)
tapply(dados$Fenois, dados$Cor, sd)
# Linear regression model- ---------------------------------------------
fit1 <- lm(Fenois ~ Cor, data = dados)
anova(fit1)
plot(residuals(fit1) ~ fitted(fit1))
plot(residuals(fit1) ~ Cor)
# Double Linear regression model ---------------------------------------
dados$id <- 1
fit2 <- mcglm(c(Fenois ~ Cor), list(mc_dglm(~ Cor, id = "id", data =
dados)),
covariance = "expm", data = dados)
summary(fit2)
summary(fit1)
summary(fit2)
cbind(coef(fit1), coef(fit2, type = "beta")$Estimates)
cbind(sqrt(diag(vcov(fit1))), coef(fit2, type = "beta", std.error =
TRUE)$Std.error)
# Example 2 ------------------------------------------------------------
dados2 <- read.table("HETE.csv", header = TRUE, sep = ";", dec = ",")
with(dados2, boxplot(y ~ x))
# Note that, the A33 has no variance, so we need to remove this level.
dados2 <- dados2[which(dados2$x != "A33"),]
dados2$x <- droplevels(dados2$x)
tapply(dados2$y, dados2$x, sd)
# Linear regression model ----------------------------------------------
fit_lm <- lm(y ~ x, data = dados2)
summary(fit_lm)
plot(residuals(fit_lm) ~ fitted(fit_lm))
plot(residuals(fit_lm) ~ dados2$x)
# Double linear regression model ---------------------------------------
dados2$id <- 1
fit_dlm <- mcglm(c(y ~ x), list(mc_dglm(~ x, id = "id", data = dados2)),
covariance = "expm", data = dados2)
summary(fit_dlm)
# Comparing estimates and standard errors ------------------------------
cbind(coef(fit_lm), coef(fit_dlm, type = "beta")$Estimates)
cbind(sqrt(diag(vcov(fit_lm))), coef(fit_dlm, type = "beta", std.error =
TRUE)$Std.error)
Em 11 de novembro de 2016 17:01, Luiz Leal <richfield1974 em yahoo.com>
escreveu:
> 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?
> Att
> Luiz
>
>
> On Friday, November 4, 2016 5:05 PM, Wagner Bonat <wbonat em gmail.com>
> wrote:
>
>
> Sim, vc pode usar qualquer dos testes padrões.
> A única diferença é que vc deve usar os beta e erros padrões que vem deste
> modelo.
> Não sei se vc pode incluir sua própria matriz de covariancia nas funções
> do R.
> Talvez, sim. Se não tem que implementar.
> O Walmes é muito bom nestes testes de comparações múltiplas talvez ele
> possa ajudar.
>
> O que houve q não conseguiu instalar o pacote?
> Tentou pelo github ou do CRAN?
>
> https://cran.r-project.org/web/packages/mcglm/index.html
> https://github.com/wbonat/mcglm
>
>
>
>
> Em 4 de novembro de 2016 19:46, Luiz Leal <richfield1974 em yahoo.com>
> escreveu:
>
> 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?
> Desde já agradeço
> Luiz
> PS. Não consegui instalar o pacote
>
>
> On Friday, November 4, 2016 4:13 PM, Wagner Bonat via R-br <
> r-br em listas.c3sl.ufpr.br> wrote:
>
>
> Caros,
>
> 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
> ser facilmente resolvido e mostrando o efeito disso no modelo.
>
> # Example 2 ------------------------------ ------------------------------
> Fenois = c(337.311, 344.874, 342.353, 325.546, 333.950, 330.588, 328.067,
> 328.067, 318.824, 331.429, 333.950, 334.790, 336.471, 338.151,
> 342.353, 259.160, 252.437, 268.403, 265.882, 266.723, 287.731,
> 88.571, 88.571, 90.252, 41.513, 52.437, 49.076, 88.571,
> 88.571, 90.252, 64.202, 60.000, 61.681)
> Cor = factor(c(rep("ambar",6), rep("ambar_claro",3), rep("ambar",6),
> rep("ambar_claro",6),rep(" branco",6),
> rep("extra_ambar_claro",3), rep("branco",3)))
>
>
> # Exploratory analysis
> boxplot(Fenois ~ Cor)
> tapply(Fenois, Cor, sd)
> dados <- data.frame(Fenois, Cor)
> dados$id <- 1
>
> # Fitting
> fit1 <- mcglm(c(Fenois ~ Cor), list(mc_id(dados)), data = dados)
> fit2 <- mcglm(c(Fenois ~ Cor), list(mc_dglm(~ Cor, id = "id", data =
> dados)),
> covariance = "expm", data = dados)
> # Goodness-of-fit
> gof(fit1)
> gof(fit2)
>
> # Comparing estimates and standard errors
> coef(fit1, type = "beta", std.error = TRUE)
> coef(fit2, type = "beta", std.error = TRUE)
>
> O interessante é que a estimativa pontual é exatamente a mesma, porém
> olha a enorme diferença nos erros padrões dos betas.
>
> --
> Wagner Hugo Bonat
> ------------------------------ ------------------------------
> ------------------------------ ----
> Department of Mathematics and Computer Science (IMADA)
> University of Southern Denmark (SDU) and
> Laboratório de Estatística e Geoinformação (LEG)
> Universidade Federal do Paraná (UFPR)
>
>
> ______________________________ _________________
> R-br mailing list
> R-br em listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/ cgi-bin/mailman/listinfo/r-br
> <https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br>
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br- guia
> <http://www.leg.ufpr.br/r-br-guia>) e forne� c�igo m�imo reproduz�el.
>
>
>
>
> --
> Wagner Hugo Bonat
> ------------------------------------------------------------
> ----------------------------------
> Department of Mathematics and Computer Science (IMADA)
> University of Southern Denmark (SDU) and
> Laboratório de Estatística e Geoinformação (LEG)
> Universidade Federal do Paraná (UFPR)
>
>
>
>
--
Wagner Hugo Bonat
----------------------------------------------------------------------------------------------
Department of Mathematics and Computer Science (IMADA)
University of Southern Denmark (SDU) and
Laboratório de Estatística e Geoinformação (LEG)
Universidade Federal do Paraná (UFPR)
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20161113/ccacbbd7/attachment-0001.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: HETE.csv
Tipo: text/csv
Tamanho: 360 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20161113/ccacbbd7/attachment-0001.csv>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: mcglm_0.4.0.tar.gz
Tipo: application/x-gzip
Tamanho: 2064531 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20161113/ccacbbd7/attachment-0001.bin>
Mais detalhes sobre a lista de discussão R-br