[R-br] componente da variância no modelo de efeitos aleatórios

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[R-br] componente da variância no modelo de efeitos aleatórios

R-br mailing list
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



_______________________________________________
R-br mailing list
[hidden email]
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.

1574880488033blob.jpg (11K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: [R-br] componente da variância no modelo de efeitos aleatórios

R-br mailing list
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.

_______________________________________________
R-br mailing list
[hidden email]
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.

1574880488033blob.jpg (11K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: [R-br] componente da variância no modelo de efeitos aleatórios

R-br mailing list
Prezado Walmes, muito obrigado.


Pelo que entendi não é possível "sigma^2_b""sigma^2_e" utilizando a função gls().

A função VarCorr(modelo2) fornece "sigma^2_b""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 <[hidden email]> 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.


_______________________________________________
R-br mailing list
[hidden email]
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.
Reply | Threaded
Open this post in threaded view
|

Re: [R-br] componente da variância no modelo de efeitos aleatórios

R-br mailing list
Respostas dentro da mensagem.

À disposição.
Walmes.

On Thu, Nov 28, 2019 at 2:49 PM Luiz Leal <[hidden email]> wrote:
Prezado Walmes, muito obrigado.


Pelo que entendi não é possível "sigma^2_b""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""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

_______________________________________________
R-br mailing list
[hidden email]
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.