Jika saya memahami pertanyaan Anda dengan benar, ini cukup mudah. Anda hanya perlu memutuskan distribusi apa yang Anda inginkan kesalahan Anda, dan menggunakan fungsi generasi acak yang sesuai.
Ada beberapa distribusi miring, jadi Anda perlu mencari tahu mana yang Anda suka. Selain itu, sebagian besar distribusi miring (misalnya, log normal, chi-squared, Gamma, Weibull, dll.) Condong ke kanan, sehingga beberapa adaptasi kecil akan diperlukan (misalnya, kalikan dengan−1).
Berikut ini contoh memodifikasi kode Anda:
set.seed(5840) # this makes the example exactly reproducible
N <- 100
x <- rnorm(N)
beta <- 0.4
errors <- rlnorm(N, meanlog=0, sdlog=1)
errors <- -1*errors # this makes them left skewed
errors <- errors - 1 # this centers the error distribution on 0
y <- 1 + x*beta + errors
Saya harus mencatat pada titik ini bahwa regresi tidak membuat asumsi tentang distribusi X atau Y, hanya tentang kesalahan, ε(lihat di sini: Bagaimana jika residu terdistribusi normal, tetapi y tidak? ). Jadi, itulah fokus jawaban saya di atas.
Pembaruan: Ini adalah versi miring kanan dengan kesalahan didistribusikan sebagai Weibull:
set.seed(5840) # this makes the example exactly reproducible
N <- 100
x <- rnorm(N)
beta <- 0.4
errors <- rweibull(N, shape=1.5, scale=1)
# errors <- -1*errors # this makes them left skewed
errors <- errors - factorial(1/1.5) # this centers the error distribution on 0
y <- 1 + x*beta + errors
Data Weibull sudah benar miring, jadi kita tidak perlu mengubah arahnya (yaitu, kita menjatuhkan -1*errors
bagian). Juga, dari halaman Wikipedia untuk distribusi Weibull, kita melihat bahwa rata-rata Weibull adalah:E[ W] = ( 1 / s h a p e ) !. Kami ingin mengurangi nilai itu dari masing-masing kesalahan sehingga distribusi kesalahan yang dihasilkan terpusat0. Itu memungkinkan bagian struktural (yaitu, 1 + x*beta
) kode Anda untuk secara akurat mencerminkan bagian struktural dari proses pembuatan data.
The distribusi ExGaussian adalah jumlah normal dan eksponensial. Ada fungsi ? RexGAUS dalam paket gamlss.dist untuk menghasilkan ini. Saya tidak memiliki paket itu, tetapi Anda harus dapat mengadaptasi kode saya di atas tanpa terlalu banyak kesulitan. Anda juga bisa menghasilkan variabel normal acak (via rnorm()
) dan eksponensial (via rexp()
) dan menjumlahkannya dengan mudah. Ingatlah untuk mengurangi mean populasi,μ + 1 / λ, dari setiap kesalahan sebelum menambahkan kesalahan ke bagian struktural dari proses pembuatan data. (Berhati-hatilah untuk tidak mengurangi mean sampelmean(errors)
!)
Beberapa komentar terakhir dan tidak terkait: Contoh kode Anda dalam pertanyaan agak kacau (artinya jangan tersinggung). Karena rnorm(N)
menghasilkan data dengan mean=0
dan sd=1
secara default, 0.4*rnorm(N)
akan menghasilkan rnorm(N, mean=0, sd=0.4)
. Kode Anda (dan mungkin pemikiran Anda) akan jauh lebih jelas jika Anda menggunakan formulasi yang terakhir. Selain itu, kode Anda untuk beta
tampaknya membingungkan. Kami biasanya memikirkanβdalam model tipe regresi sebagai parameter, bukan variabel acak. Artinya, ini adalah konstanta yang tidak diketahui yang mengatur perilaku proses pembuatan data, tetapi sifat stokastik dari proses tersebut dienkapsulasi oleh kesalahan. Ini bukan cara kami memikirkannya ketika kami bekerja dengan model multilevel, dan kode Anda tampaknya setengah antara model regresi standar dan kode untuk model regresi multilevel. Menentukan beta Anda secara terpisah adalah ide yang bagus untuk menjaga kejelasan konseptual kode, tetapi untuk model regresi standar, Anda hanya akan menetapkan satu angka untuk setiap beta (misalnya, beta0 <- 1; beta1 <- .04
).