componente da variância no modelo de efeitos aleatórios
Prezados, boa tarde. No modelo de efeitos aleatórios estimamos a componente da variância conforme abaixo. y <-c(2370, 1687, 2592, 2283, 2910, 3020, 1282, 1527, 871, 1025, 825, 920, 562, 321, 636, 317, 485, 842, 173, 127, 132, 150, 129, 227, 193, 71, 82, 62, 96, 44)x <- as.factor(c("A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D", "E", "E", "E", "E", "E", "E"))summary(aov(y~x)) Mas não sei como calcular quando utilizo mínimos quadrados generalizados. y <-c(2370, 1687, 2592, 2283, 2910, 3020, 1282, 1527, 871, 1025, 825, 920, 562, 321, 636, 317, 485, 842, 173, 127, 132, 150, 129, 227, 193, 71, 82, 62, 96, 44) x <- as.factor(c("A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D", "E", "E", "E", "E", "E", "E"))library(nlme) modelo <- gls(y~x, weights = varExp() , method = "ML") Alguém sabe como calculo as componentes da variância neste caso? Desde já agradeço Luiz
Luiz, Você terá que mudar a função que está usando. library(nlme) table(x) lattice::xyplot(y ~ x) modelo0 <- aov(y ~ x) a <- anova(modelo0) a # Componentes de variância pelo método ANOVA. c("sigma^2_b" = (a["x", "Mean Sq"] - a["Residuals", "Mean Sq"])/6, "sigma^2_e" = a["Residuals", "Mean Sq"]) modelo1 <- lme(y ~ 1, random = ~1 | x, method = "REML") modelo1 # Componentes de variância pelo método REML. VarCorr(modelo1) modelo2 <- lme(y ~ 1, random = ~1 | x, weights = varExp(), method = "ML") modelo2 # Componentes de variância pelo método REML. # ATTENTION: não são imediatamente comparáveis com os anteriores. VarCorr(modelo2) À disposição. Walmes.
Prezado Walmes, muito obrigado. Pelo que entendi não é possível "sigma^2_b" e "sigma^2_e" utilizando a função gls(). A função VarCorr(modelo2) fornece "sigma^2_b" e "sigma^2_e", correto? Variance StdDev (Intercept) 667641.149 817.09311 Residual 6041.696 77.72834 Tenho uma dúvida: a variância do intercepto é "sigma^2_b"? Muito obrigado. Sua contribuição está sendo muito válida no meu trabalho.Luiz On Wednesday, November 27, 2019, 07:13:37 PM GMT-2, Walmes Zeviani <walmeszeviani@gmail.com> wrote: Luiz, Você terá que mudar a função que está usando. library(nlme) table(x) lattice::xyplot(y ~ x) modelo0 <- aov(y ~ x) a <- anova(modelo0) a # Componentes de variância pelo método ANOVA. c("sigma^2_b" = (a["x", "Mean Sq"] - a["Residuals", "Mean Sq"])/6, "sigma^2_e" = a["Residuals", "Mean Sq"]) modelo1 <- lme(y ~ 1, random = ~1 | x, method = "REML") modelo1 # Componentes de variância pelo método REML. VarCorr(modelo1) modelo2 <- lme(y ~ 1, random = ~1 | x, weights = varExp(), method = "ML") modelo2 # Componentes de variância pelo método REML. # ATTENTION: não são imediatamente comparáveis com os anteriores. VarCorr(modelo2) À disposição.Walmes.
Respostas dentro da mensagem. À disposição. Walmes. On Thu, Nov 28, 2019 at 2:49 PM Luiz Leal <richfield1974@yahoo.com> wrote:
Prezado Walmes, muito obrigado.
Pelo que entendi não é possível "sigma^2_b" e "sigma^2_e" utilizando a função gls().
Não. Tem-se que usar um modelo de efeitos aleatórios.
A função VarCorr(modelo2) fornece "sigma^2_b" e "sigma^2_e", correto?
Sim. Mas a interpretação é mais específica. "sigma^2_e" é a variância do erro ao redor do ponto médio da curva em um ponto específico do suporte da covariável que eu imagino ser o 0. A variância em outros pontos da curva será outra, obviamente, porque você declarou um modelo heterocedástico. # Média ajustada. m <- unique(fitted(modelo2)) # Modelo ajustado. plot(y ~ as.integer(x)) lines(m ~ seq_along(m)) # Variância estimada é função da média. # s2(v) = exp(2* t * v) var_estim <- 6041.676 * exp(2 * 0.0009743905 * m) var_amost <- tapply(y, x, var) cbind(var_amost, var_estim) # Resíduos crus e padronizados. plot(residuals(modelo2) ~ fitted(modelo2)) plot(residuals(modelo2, type = "pearson") ~ fitted(modelo2)) # Calculando o resíduo padronizado. r_my <- residuals(modelo2)/ sqrt(6041.676 * exp(2 * 0.0009743905 * fitted(modelo2))) # Comparação. cbind(r_my, residuals(modelo2, type = "pearson"))
Variance StdDev (Intercept) 667641.149 817.09311 Residual 6041.696 77.72834
Tenho uma dúvida: a variância do intercepto é "sigma^2_b"?
Sim.
Muito obrigado. Sua contribuição está sendo muito válida no meu trabalho. Luiz
Muito obrigado Walmes. On Monday, December 2, 2019, 12:40:19 PM GMT-2, Walmes Zeviani <walmeszeviani@gmail.com> wrote: Respostas dentro da mensagem. À disposição.Walmes. On Thu, Nov 28, 2019 at 2:49 PM Luiz Leal <richfield1974@yahoo.com> wrote: Prezado Walmes, muito obrigado. Pelo que entendi não é possível "sigma^2_b" e "sigma^2_e" utilizando a função gls(). Não. Tem-se que usar um modelo de efeitos aleatórios. A função VarCorr(modelo2) fornece "sigma^2_b" e "sigma^2_e", correto? Sim. Mas a interpretação é mais específica. "sigma^2_e" é a variância do erro ao redor do ponto médio da curva em um ponto específico do suporte da covariável que eu imagino ser o 0. A variância em outros pontos da curva será outra, obviamente, porque você declarou um modelo heterocedástico. # Média ajustada. m <- unique(fitted(modelo2)) # Modelo ajustado. plot(y ~ as.integer(x)) lines(m ~ seq_along(m)) # Variância estimada é função da média. # s2(v) = exp(2* t * v) var_estim <- 6041.676 * exp(2 * 0.0009743905 * m) var_amost <- tapply(y, x, var) cbind(var_amost, var_estim) # Resíduos crus e padronizados. plot(residuals(modelo2) ~ fitted(modelo2)) plot(residuals(modelo2, type = "pearson") ~ fitted(modelo2)) # Calculando o resíduo padronizado. r_my <- residuals(modelo2)/ sqrt(6041.676 * exp(2 * 0.0009743905 * fitted(modelo2))) # Comparação. cbind(r_my, residuals(modelo2, type = "pearson")) Variance StdDev (Intercept) 667641.149 817.09311 Residual 6041.696 77.72834 Tenho uma dúvida: a variância do intercepto é "sigma^2_b"? Sim. Muito obrigado. Sua contribuição está sendo muito válida no meu trabalho.Luiz
participantes (2)
-
Luiz Leal -
Walmes Zeviani