Simule regressão linear com heterocedasticidade

9

Estou tentando simular um conjunto de dados que corresponda aos dados empíricos que tenho, mas não tenho certeza de como estimar os erros nos dados originais. Os dados empíricos incluem heterocedasticidade, mas não estou interessado em transformá-los, mas sim usando um modelo linear com um termo de erro para reproduzir simulações dos dados empíricos.

Por exemplo, digamos que eu tenho um conjunto de dados empíricos e um modelo:

n=rep(1:100,2)
a=0
b = 1
sigma2 = n^1.3
eps = rnorm(n,mean=0,sd=sqrt(sigma2))
y=a+b*n + eps
mod <- lm(y ~ n)

usando plot(n,y), obtemos o seguinte. insira a descrição da imagem aqui

No entanto, se eu tentar simular os dados simulate(mod), a heterocedasticidade é removida e não capturada pelo modelo.

Eu posso usar um modelo de mínimos quadrados generalizado

VMat <- varFixed(~n)
mod2 = gls(y ~ n, weights = VMat)

que fornece um melhor ajuste do modelo baseado no AIC, mas não sei como simular dados usando a saída.

Minha pergunta é: como crio um modelo que me permita simular dados para corresponder aos dados empíricos originais (n e y acima). Especificamente, preciso de uma maneira de estimar o sigma2, o erro, usando um modelo?

user44796
fonte
11
Portanto, o modelo linear não captura a heterocedasticidade condicional, a menos que tente explicitamente fazê-lo, usando uma das poucas abordagens. As técnicas econométricas padrão ajustam os erros padrão nos parâmetros para explicar a heterocedasticidade, mas não o modelam explicitamente.
generic_user
Você está certo. Estou tentando usar um modelo linear para capturar a heterogeneidade. Eu acho que eu deveria estar usando um modelo de mínimos quadrados generalizado. Se houver outras recomendações, tentarei.
user44796
Existe um erro no seu código, você deve usar `lm (y ~ n)`
kjetil b halvorsen
11
Eu não entendo sua pergunta, porque seu código realiza exatamente o que você parece estar pedindo no título: simula uma regressão linear com erros heterocedásticos. Você está pedindo métodos para estimar algum tipo de modelo para a heterocedasticidade? Nesse caso, você precisa especificar um modelo!
whuber
Espero ter esclarecido minha pergunta com edições. Na pergunta acima, n e y representam os dados empíricos. Desejo ajustar um modelo aos dados e depois usá-lo para gerar dados simulados que correspondam à média e aos resíduos dos dados originais.
user44796

Respostas:

9

Para simular dados com uma variação de erro variável, é necessário especificar o processo de geração de dados para a variação de erro. Como foi apontado nos comentários, você fez isso quando gerou seus dados originais. Se você possui dados reais e deseja tentar isso, basta identificar a função que especifica como a variação residual depende de suas covariáveis. A maneira padrão de fazer isso é ajustar o seu modelo, verificar se é razoável (além da heterocedasticidade) e salvar os resíduos. Esses resíduos se tornam a variável Y de um novo modelo. Abaixo, fiz isso para o seu processo de geração de dados. (Não vejo onde você define a semente aleatória, portanto, esses não serão literalmente os mesmos dados, mas devem ser semelhantes, e você pode reproduzir os meus exatamente usando minha semente.)

set.seed(568)  # this makes the example exactly reproducible

n      = rep(1:100,2)
a      = 0
b      = 1
sigma2 = n^1.3
eps    = rnorm(n,mean=0,sd=sqrt(sigma2))
y      = a+b*n + eps
mod    = lm(y ~ n)
res    = residuals(mod)

windows()
  layout(matrix(1:2, nrow=2))
  plot(n,y)
  abline(coef(mod), col="red")
  plot(mod, which=3)

insira a descrição da imagem aqui

Observe que Rs ? Plot.lm fornecerá um gráfico (cf., aqui ) da raiz quadrada dos valores absolutos dos resíduos, sobrepostos de maneira útil com um ajuste inferior, que é exatamente o que você precisa. (Se você tiver várias covariáveis, poderá avaliar isso em relação a cada covariada separadamente.) Há a menor sugestão de curva, mas parece que uma linha reta faz um bom trabalho ao ajustar os dados. Então, vamos ajustar explicitamente esse modelo:

res.mod = lm(sqrt(abs(res))~fitted(mod))
summary(res.mod)
# Call:
# lm(formula = sqrt(abs(res)) ~ fitted(mod))
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -3.3912 -0.7640  0.0794  0.8764  3.2726 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 1.669571   0.181361   9.206  < 2e-16 ***
# fitted(mod) 0.023558   0.003157   7.461 2.64e-12 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.285 on 198 degrees of freedom
# Multiple R-squared:  0.2195,  Adjusted R-squared:  0.2155 
# F-statistic: 55.67 on 1 and 198 DF,  p-value: 2.641e-12
windows()
  layout(matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
  plot(res.mod, which=1)
  plot(res.mod, which=2)
  plot(res.mod, which=3)
  plot(res.mod, which=5)

insira a descrição da imagem aqui

Não precisamos nos preocupar que a variação residual pareça estar aumentando também no gráfico de localização em escala para esse modelo - isso essencialmente tem que acontecer. Há novamente a menor sugestão de curva, para que possamos tentar ajustar um termo ao quadrado e ver se isso ajuda (mas não ajuda):

res.mod2 = lm(sqrt(abs(res))~poly(fitted(mod), 2))
summary(res.mod2)
# output omitted
anova(res.mod, res.mod2)
# Analysis of Variance Table
# 
# Model 1: sqrt(abs(res)) ~ fitted(mod)
# Model 2: sqrt(abs(res)) ~ poly(fitted(mod), 2)
#   Res.Df    RSS Df Sum of Sq     F Pr(>F)
# 1    198 326.87                          
# 2    197 326.85  1  0.011564 0.007 0.9336

Se estivermos satisfeitos com isso, agora podemos usar esse processo como um complemento para simular dados.

set.seed(4396)  # this makes the example exactly reproducible
x = n
expected.y = coef(mod)[1] + coef(mod)[2]*x
sim.errors = rnorm(length(x), mean=0,
                   sd=(coef(res.mod)[1] + coef(res.mod)[2]*expected.y)^2)
observed.y = expected.y + sim.errors

Observe que esse processo não tem mais garantia de encontrar o verdadeiro processo de geração de dados do que qualquer outro método estatístico. Você usou uma função não linear para gerar os SDs de erro e nós a aproximamos com uma função linear. Se você realmente conhece o verdadeiro processo de geração de dados a priori (como neste caso, porque simulou os dados originais), você também pode usá-lo. Você pode decidir se a aproximação aqui é boa o suficiente para seus propósitos. Entretanto, normalmente não conhecemos o verdadeiro processo de geração de dados e, com base no barbeador da Occam, executamos a função mais simples que se encaixa adequadamente aos dados que fornecemos à quantidade de informações disponíveis. Você também pode tentar splines ou abordagens mais sofisticadas, se preferir. As distribuições bivariadas são razoavelmente semelhantes a mim,

insira a descrição da imagem aqui

- Reinstate Monica
fonte
Esta era realmente uma conclusão que eu estava começando a chegar, mas nunca teria chegado a uma resposta tão elegante.
user44796
5

Você precisa modelar a heterocedasticidade. Uma abordagem é através do pacote R (CRAN) dglm, modelo linear generalizado de dispersão. Esta é uma extensão de glm's que, além do habitual glm, encaixa um segundo glm para dispersão dos resíduos do primeiro glm. Não tenho experiência com esses modelos, mas eles parecem promissores ... Aqui estão alguns códigos:

n <- rep(1:100,2)
a <- 0
b <- 1
sigma2 <- n^1.3
eps <- rnorm(n,mean=0,sd=sqrt(sigma2))
y <- a+b*n + eps
mod <- lm(y ~ n)

library(dglm)  ### double glm's

mod2   <-  dglm(y ~ n, ~ n, gaussian,ykeep=TRUE,xkeep=TRUE,zkeep=TRUE)
### This uses log link for the dispersion part, should also try identity link ..

y2 <-  simulate(mod2)

plot(n, y2$sim_1)

mod3  <-  dglm(y ~ n, ~ n, gaussian, dlink="identity", ykeep=TRUE,xkeep=TRUE,zkeep=TRUE)  ### This do not work because it leads to negative weights!

O gráfico simulado é mostrado abaixo:

insira a descrição da imagem aqui

O gráfico parece que a simulação usou a variação estimada, mas não tenho certeza, pois a função simulate () não possui métodos para os dglm ...

(Outra possibilidade de examinar é usar o Rpacote gamlss, que usa outra abordagem para modelar a variação em função das covariáveis.)

kjetil b halvorsen
fonte
11
o modelo linear generalizado duplo parece modelar os dados originais adequadamente. Não estou claro como o erro residual é modelado usando o predizer (). Vou ter que dar uma olhada nisso.
user44796