3 votos

Generación de conjuntos de datos aleatorios para la regresión lineal con pendiente aleatoria y término de error en R

Quiero probar los efectos del tamaño de la muestra en la regresión de Deming utilizando datos emparejados simulados en R.

Como los datos están emparejados, el valor esperado de la pendiente debe ser 1 y el intercepto 0.

El código que tengo hasta ahora es:

x <- runif(16,0,25)
y <- rnorm(16,1,0.05)*x + rnorm(16,0,0.5)

Esto genera 16 pares de datos aleatorios entre 0 y 25 con pendiente aleatoria sobre 1 y error aleatorio sobre 0 que es lo que quiero.

Mis preguntas son:

1) ¿Qué desviaciones estándar debo elegir para que la pendiente y el término de error caigan para dar un alfa de 0,05 (95% de confianza). ¿Debería ser 1,96 para el término de error?

2) ¿Debo incluir también un intercepto aleatorio?

1voto

TheEvilPenguin Puntos 3904

Ahora mismo tu código está generando 16 puntos de datos correlacionados, aunque tu terminología para algunos de los que estás hablando no es del todo correcta.

No estás generando una "pendiente aleatoria" en el sentido de lo que la gente habla cuando se refiere a los modelos mixtos (por ejemplo, interceptos aleatorios y pendientes aleatorias que varían por alguna variable de agrupación). En este momento está generando un único conjunto de datos con una única pendiente, pero en el que ha extraído la pendiente aleatoriamente de una distribución normal con una media de 1 y una desviación estándar de 0,05. Es decir, su pendiente será cercana a 1, pero será un poco diferente cada vez que vuelva a ejecutar esas líneas de código, a menos que vuelva a establecer la semilla aleatoria al mismo valor antes de cada ejecución del código (por ejemplo, set.seed(101) antes de cada vez que lo ejecute -- entonces obtendrá el mismo valor de pendiente cada vez).

Además, no estás dibujando un "error aleatorio" de alrededor de cero. El error que estás añadiendo a tu variable y se extrae aleatoriamente de una distribución normal, con una media de cero y una desviación estándar de 0,5 (no es exactamente cero, pero es menor que tu pendiente).

Además, no estoy seguro de lo que quieres decir al especificar las desviaciones estándar para obtener un alfa de 0,05. Un nivel alfa es una propiedad de una prueba estadística de hipótesis nula, no de un conjunto de datos.

Lo que puede hacer en esta situación es estimar el intervalo de confianza del 95% para su pendiente. Esto le dirá si su pendiente es, con un 95% de confianza, diferente de cero y le dará un intervalo sobre el que puede utilizar para estimar los valores probables de la pendiente verdadera (parámetro poblacional), dado el conjunto de datos de la muestra. El IC será una función del error estimado de la pendiente y del tamaño de la muestra. El error estándar de una estimación de la pendiente viene dado por:

$se\hat{\beta}_1 = \sqrt{\frac{\frac{1}{n-2} \sum_{i=1}^n (\hat{y}_i-y_i)^2} {\sum_{i=1}^n (x_i - \bar{x}_i)^2}}$

A continuación, puede convertir esto en un intervalo de confianza para la pendiente como $b_i \pm t_{\alpha/2, n-2}\times se\hat{\beta}_1 $ , donde $b_i$ es la pendiente que has estimado a partir de tu modelo, y $t_{\alpha/2, n-2}$ es el valor t crítico para su alfa especificado, con n-2 grados de libertad.

Así, con una sola muestra como ésta, se puede estimar el intercepto y su error, y la pendiente de la muestra y su error. Cada uno de ellos te dará intervalos de confianza para decirte si es probable que el verdadero intercepto y la pendiente de la población que estás tratando de estimar sean distintos de cero. Recuerde que la potencia estadística disminuye con el tamaño de la muestra, por lo que con muestras pequeñas como la que ha elegido, los intervalos de confianza probablemente incluyan el cero, especialmente si su error es del orden de la pendiente.

Si está tratando de probar la potencia y/o la tasa de error de tipo I, puede generar múltiples muestras para varios ajustes de la pendiente, la varianza del error y el tamaño de la muestra. Puede probar este código:

    # Set your random generator seed for reproducibility
    set.seed(101)

    # Some parameters to try
    n_obs <- c(10, 15, 20, 25)
    slopes <- c(0, 1, 2, 3)
    stdevs <- c(.5, 1, 5, 10)
    intercept <- 0 # Constant for all sims

    # Generate a df with all combinations of the parameters
    parameter_matrix <- expand.grid(N = n_obs, 
                                    Intercept = intercept, 
                                    Slope = slopes, 
                                    SD = stdevs)

    # Create some columns in the matrix to store model outputs
    parameter_matrix$Intercept_estimate <- NA
    parameter_matrix$Intercept_SE <- NA
    parameter_matrix$Intercept_t <- NA
    parameter_matrix$Intercept_p <- NA
    parameter_matrix$Slope_estimate <- NA
    parameter_matrix$Slope_SE <- NA
    parameter_matrix$Slope_t <-NA
    parameter_matrix$Slope_p <- NA

    # Loop over each combination of parameters and generate a sample, fit a model, 
    # and save the estimates back to parameter_matrix
    for (i in 1:nrow(parameter_matrix)){

      # Generate data
      x <- runif(parameter_matrix$N[i], 0, 20)
      y <- x*rnorm(parameter_matrix$N[i], 
                   mean = parameter_matrix$Slope[i], 
                   sd = parameter_matrix$SD[i])

      # Fit a model to the data
      mod <- summary(lm(y ~ x))

      # Pull out the parameter estimates, errors, etc, and save them
      parameter_matrix$Intercept_estimate[i] <- mod$coefficients[1,1]
      parameter_matrix$Intercept_SE[i] <- mod$coefficients[1, 2]
      parameter_matrix$Intercept_t[i] <- mod$coefficients[1, 3]
      parameter_matrix$Intercept_p[i] <- mod$coefficients[1, 4]
      parameter_matrix$Slope_estimate[i] <- mod$coefficients[2, 1]
      parameter_matrix$Slope_SE[i] <- mod$coefficients[2, 2]
      parameter_matrix$Slope_t[i] <- mod$coefficients[2, 3]
      parameter_matrix$Slope_p[i] <- mod$coefficients[2, 4]
    }

A partir de aquí puede ver cómo sus estimaciones de la $\hat{\beta}$ y su cambio de error en función de la pendiente especificada, la varianza del error y el tamaño de la muestra (por ejemplo, con gráficos de dispersión que muestren comparaciones por pares de los parámetros, con otros parámetros delineados por color y/o forma). Para ser aún más completo, querrá incrustar esto en otro bucle para generar múltiples (cientos o miles) de simulaciones para cada ajuste de parámetros para tener una idea de la variabilidad del muestreo aleatorio. Recuerde que ocasionalmente verá grandes valores t/pequeños valores p, incluso con una pendiente de 0 (es decir, errores de tipo I), porque esa es la naturaleza del muestreo. También verá pendientes no significativas, incluso cuando especifique una pendiente distinta de cero (error de tipo II), especialmente con muestras pequeñas y/o una mayor varianza de error.

Las pendientes aleatorias (que es un punto sobre el que preguntas en (2)) no son apropiadas para muestras extraídas de una única población como ésta. Las pendientes aleatorias son apropiadas cuando la muestra tiene una estructura anidada, como si se trata de un muestreo de estudiantes dentro de las aulas. Por ejemplo, usted examina a 20 estudiantes en cada una de las 20 aulas, por lo que tiene 20 x 20 = 400 muestras, pero esas 400 muestras tienen una estructura correlacional anidada (se espera que los estudiantes dentro de un aula sean más similares entre sí que los estudiantes extraídos de diferentes aulas).

Se querrá ajustar un intercepto aleatorio para cada aula, además de estimar el intercepto fijo (promedio) entre las aulas, más las pendientes para cualquier otro efecto en el que esté interesado (por ejemplo, el efecto del estatus socioeconómico (variable x) en alguna puntuación de la prueba (variable y)). En este caso, debería considerar también el ajuste de las pendientes aleatorias para cada variable de agrupación (aula), dependiendo de su diseño.

Información sobre cómo simular datos para un diseño de efectos mixtos se puede encontrar aquí.

Es un poco más complicado que generar dos vectores de datos correlacionados.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X