Mensimulasikan regresi linier dengan heteroskedastisitas

9

Saya mencoba mensimulasikan dataset yang cocok dengan data empiris yang saya miliki, tetapi saya tidak yakin bagaimana memperkirakan kesalahan dalam data asli. Data empiris mencakup heteroskedastisitas, tetapi saya tidak tertarik untuk mengubahnya, tetapi menggunakan model linier dengan istilah kesalahan untuk mereproduksi simulasi data empiris.

Sebagai contoh, katakanlah saya memiliki beberapa dataset empiris dan model:

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)

menggunakan plot(n,y)kita mendapatkan yang berikut ini. masukkan deskripsi gambar di sini

Namun, jika saya mencoba mensimulasikan data, simulate(mod)heteroskedastisitas dihapus dan tidak ditangkap oleh model.

Saya dapat menggunakan model kuadrat terkecil umum

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

yang memberikan model fit yang lebih baik berdasarkan AIC, tapi saya tidak tahu bagaimana mensimulasikan data menggunakan output.

Pertanyaan saya adalah, bagaimana cara membuat model yang akan memungkinkan saya untuk mensimulasikan data agar sesuai dengan data asli, empiris (n dan y di atas). Secara khusus, saya butuh cara untuk memperkirakan sigma2, kesalahan, menggunakan salah satu menggunakan model?

pengguna44796
sumber
1
Jadi model linier tidak akan menangkap heteroskedastisitas bersyarat kecuali ia secara eksplisit mencoba melakukannya, menggunakan salah satu dari beberapa pendekatan. Teknik ekonometrik standar menyesuaikan kesalahan standar pada parameter untuk memperhitungkan heteroskedastisitas, tetapi mereka tidak secara eksplisit memodelkannya.
generic_user
Kamu benar. Saya mencoba menggunakan model linier untuk menangkap heterogenitas. Saya pikir saya harus menggunakan model kuadrat terkecil umum. Jika ada rekomendasi lain, saya akan mencobanya.
user44796
Tidak Ada KESALAHAN DALAM KODE ANDA, ANDA HARUS MENGGUNAKAN `lm (y ~ n)`
kjetil b halvorsen
1
Saya tidak mengerti pertanyaan Anda, karena kode Anda memenuhi apa yang tampaknya Anda tanyakan dalam judulnya: ini mensimulasikan regresi linier dengan kesalahan heteroskedastik. Apakah Anda meminta metode untuk memperkirakan semacam model untuk heteroskedastisitas? Jika demikian, maka Anda perlu menentukan model!
whuber
Semoga saya telah mengklarifikasi pertanyaan saya dengan suntingan. Dalam pertanyaan di atas, n dan y mewakili data empiris. Saya ingin menyesuaikan model dengan data dan kemudian menggunakan model untuk menghasilkan data simulasi yang cocok dengan rata-rata dan residu dari data asli.
user44796

Jawaban:

9

Untuk mensimulasikan data dengan varians kesalahan yang bervariasi, Anda perlu menentukan proses pembuatan data untuk varians kesalahan. Seperti yang telah ditunjukkan dalam komentar, Anda melakukan itu ketika Anda menghasilkan data asli Anda. Jika Anda memiliki data nyata dan ingin mencoba ini, Anda hanya perlu mengidentifikasi fungsi yang menentukan bagaimana varians residual tergantung pada kovariat Anda. Cara standar untuk melakukan itu adalah agar sesuai dengan model Anda, periksa apakah itu masuk akal (selain heteroskedastisitas), dan simpan residu. Residu tersebut menjadi variabel Y dari model baru. Di bawah ini saya telah melakukannya untuk proses pembuatan data Anda. (Saya tidak melihat di mana Anda mengatur benih acak, jadi ini tidak akan benar-benar menjadi data yang sama, tetapi harus serupa, dan Anda dapat mereproduksi milik saya dengan menggunakan benih saya.)

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)

masukkan deskripsi gambar di sini

Perhatikan bahwa R' plot.lm' akan memberi Anda plot (lih., Di sini ) dari akar kuadrat dari nilai absolut residu, yang dilapis dengan bantuan lowess, yang tepat seperti yang Anda butuhkan. (Jika Anda memiliki banyak kovariat, Anda mungkin ingin menilai ini terhadap masing-masing kovariat secara terpisah.) Ada sedikit petunjuk kurva, tetapi ini terlihat seperti garis lurus yang berfungsi baik dalam menyesuaikan data. Jadi mari kita secara eksplisit menyesuaikan model itu:

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)

masukkan deskripsi gambar di sini

Kita tidak perlu khawatir bahwa varians residual tampaknya meningkat dalam plot skala lokasi untuk model ini juga — yang pada dasarnya harus terjadi. Ada lagi sedikit tanda kurva, sehingga kita dapat mencoba menyesuaikan istilah kuadrat dan melihat apakah itu membantu (tetapi tidak):

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

Jika kami puas dengan ini, kami sekarang dapat menggunakan proses ini sebagai tambahan untuk mensimulasikan data.

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

Perhatikan bahwa proses ini tidak lagi dijamin untuk menemukan proses pembuatan data yang sebenarnya daripada metode statistik lainnya. Anda menggunakan fungsi non-linear untuk menghasilkan SD kesalahan, dan kami memperkirakannya dengan fungsi linear. Jika Anda benar-benar mengetahui proses pembuatan data sebenarnya a-priori (seperti dalam kasus ini, karena Anda mensimulasikan data asli), Anda sebaiknya menggunakannya. Anda dapat memutuskan apakah perkiraan di sini cukup baik untuk tujuan Anda. Namun, kami biasanya tidak mengetahui proses pembuatan data yang sebenarnya, dan berdasarkan pada pisau Occam, gunakan fungsi paling sederhana yang cukup sesuai dengan data yang kami berikan pada jumlah informasi yang tersedia. Anda juga dapat mencoba pendekatan splines atau pelamun jika Anda mau. Distribusi bivariat terlihat cukup mirip dengan saya,

masukkan deskripsi gambar di sini

gung - Pasang kembali Monica
sumber
Ini sebenarnya kesimpulan yang saya mulai tuju, tetapi tidak akan pernah sampai pada jawaban yang begitu elegan.
user44796
5

Anda perlu memodelkan heteroskedastisitas. Salah satu pendekatan adalah melalui paket R (CRAN) dglm, model linear umum dispersi. Ini adalah perpanjangan dari glm yang, di samping biasa glm, cocok glm kedua untuk dispersi dari residu dari glm pertama. Saya tidak punya pengalaman dengan model seperti itu, tetapi mereka tampak menjanjikan ... Berikut adalah beberapa kode:

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!

Plot simulasi ditunjukkan di bawah ini:

masukkan deskripsi gambar di sini

Plot memang terlihat seperti simulasi telah menggunakan varians yang diperkirakan, tapi saya tidak yakin, karena fungsi mensimulasikan () tidak memiliki metode untuk dglm ...

(Kemungkinan lain untuk melihat ke dalam, adalah menggunakan Rpaket gamlss, yang menggunakan pendekatan lain untuk memodelkan varians sebagai fungsi dari kovariabel.)

kjetil b halvorsen
sumber
1
model linier ganda digeneralisasi untuk memodelkan data asli secara memadai. Saya tidak jelas bagaimana kesalahan residual dimodelkan menggunakan predict (). Saya harus melihat itu.
user44796