Si entiendo su pregunta correctamente, esto es bastante fácil. Solo necesita decidir qué distribución desea que tengan sus errores y usar la función de generación aleatoria correspondiente.
Hay varias distribuciones sesgadas, por lo que debe averiguar cuál le gusta. Además, la mayoría de las distribuciones sesgadas (p. Ej., Log normal, chi-cuadrado, gamma, Weibull, etc.) están sesgadas a la derecha, por lo que serían necesarias algunas adaptaciones menores (p. Ej., Multiplicar por ). - 1
Aquí hay un ejemplo que modifica su código:
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
Debo señalar en este punto que la regresión no hace suposiciones sobre las distribuciones de o , solo sobre los errores, (ver aquí: ¿Qué pasa si los residuos están normalmente distribuidos, pero y no lo es? ). Por lo tanto, ese fue el enfoque de mi respuesta anterior. XYε
Actualización: Aquí hay una versión sesgada a la derecha con los errores distribuidos como 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
Los datos de Weibull ya están sesgados, por lo que no necesitamos cambiar su dirección (es decir, descartamos la -1*errors
parte). Además, desde la página de Wikipedia para la distribución de Weibull, vemos que la media de un Weibull debería ser:. Queremos restar ese valor de cada uno de los errores para que la distribución del error resultante se centre en . Eso permite que la parte estructural (es decir, ) de su código refleje con precisión la parte estructural del proceso de generación de datos. mi[ W] = ( 1 / s h a p e ) !0 01 + x*beta
La distribución exgaussiana es la suma de una normal y una exponencial. Hay una función ? RexGAUS en el paquete gamlss.dist para generarlos. No tengo ese paquete, pero deberías poder adaptar mi código anterior sin demasiada dificultad. También podría generar una variable normal aleatoria (vía rnorm()
) y una exponencial (vía rexp()
) y sumarlas con bastante facilidad. Solo recuerde restar la media de la población, , de cada error antes de agregar los errores a la parte estructural del proceso de generación de datos. (¡Pero tenga cuidado de no restar la media muestral !) μ + 1 / λmean(errors)
Algunos comentarios finales no relacionados: su código de ejemplo en la pregunta está algo confuso (es decir, sin ofender). Porque rnorm(N)
genera datos con mean=0
y sd=1
por defecto, 0.4*rnorm(N)
generará rnorm(N, mean=0, sd=0.4)
. Su código (y posiblemente su pensamiento) será mucho más claro si usa la última formulación. Además, su código para beta
parece confundido. Generalmente pensamos en elβen un modelo de regresión como parámetro, no como una variable aleatoria. Es decir, es una constante desconocida que gobierna el comportamiento del proceso de generación de datos, pero la naturaleza estocástica del proceso está encapsulada por los errores. Esta no es la forma en que pensamos cuando trabajamos con modelos multinivel, y su código parece estar a medio camino entre un modelo de regresión estándar y el código para un modelo de regresión multinivel. Especificar sus betas por separado es una buena idea para mantener la claridad conceptual del código, pero para un modelo de regresión estándar, simplemente asignaría un número único a cada beta (por ejemplo, beta0 <- 1; beta1 <- .04
).