Pregunta uno

Se va a iniciar construyendo las bases de datos requeridas para este punto.

# Base de datos A
x1 <- c('A', 'B', 'C')
x2 <- 1:3
bdA <- data.frame(x1, x2)
bdA

##   x1 x2
## 1  A  1
## 2  B  2
## 3  C  3

# Base de datos B
x1 <- c('A', 'B', 'D')
x3 <- c(TRUE, FALSE, TRUE)
bdB <- data.frame(x1, x3)
bdB

##   x1    x3
## 1  A  TRUE
## 2  B FALSE
## 3  D  TRUE

Pregunta dos

Se va a proceder con la construcción de las bases de datos que se quieren.

• Complemento de b - a. Esta base de datos va a tener todos los elementos de la base de datos a y solo toma de b aquellos elementos en común con la base de datos a según la columna de enlace. Si en a hay alguna fila que no tiene elemento común con la columna de enlace con b, entonces la característica que se trae de b se queda como un NA.

complemento_b_a <- left_join(x = bdA, y = bdB, by = 'x1')
complemento_b_a

##   x1 x2    x3
## 1  A  1  TRUE
## 2  B  2 FALSE
## 3  C  3    NA

• Complemento de a - b. Esta base de datos va a tener todos los elementos de la base de datos b y solo toma de a aquellos elementos en común con la base de datos b según la columna de enlace. Si en b hay alguna fila que no tiene elemento común con la columna de enlace con a, entonces la característica que se trae de a se queda como un NA.

complemento_a_b <- right_join(x = bdA, y = bdB, by = 'x1')
complemento_a_b

##   x1 x2    x3
## 1  A  1  TRUE
## 2  B  2 FALSE
## 3  D NA  TRUE

Nótese que con este resultado, si bien toma como marco de datos a referencia a b, conserva el orden de las columnas según como fueron introducidas. Vale la pena ver si el resultado es igual con un left_join() modificado:

verificacion <- left_join(x = bdB, y = bdA, by = 'x1')
verificacion

##   x1    x3 x2
## 1  A  TRUE  1
## 2  B FALSE  2
## 3  D  TRUE NA

Nótese que se obtienen las mismas columnas, solo que en orden invertido y tal como se desea según la imagen dada.

• Intersección de a y b. Solo se toman aquellos elementos para los cuales se tienen datos completos para todas las observaciones comunes para ambas bases de datos de acuerdo con la columna de enlace.

interseccion <- inner_join(x = bdA, y = bdB, by  = 'x1')
interseccion

##   x1 x2    x3
## 1  A  1  TRUE
## 2  B  2 FALSE

• Unión. Se toman todos los elementos de las dos bases de datos y se le asignan a los valores faltantes NA.

union <- full_join(x = bdA, y = bdB, by = 'x1')
union

##   x1 x2    x3
## 1  A  1  TRUE
## 2  B  2 FALSE
## 3  C  3    NA
## 4  D NA  TRUE

Pregunta tres

En el contexto de la programación, una cheat sheet, que se traduciría literalmente como hoja de engaño es lo que se conoce coloquialmente como pastel en Antioquia o chuleta en España. Son documentos de una o dos hojas por lo general con claves de algún paquete o módulo que explica de forma sucinta como usar algunas de las funciones o métodos más importantes.

Pregunta cuatro

El cheat sheet del paquete dplyr puede ser encontrado en el siguiente enlace: clic aquí

Pregunta cinco

Al nombrar cada una de las posibles bases de datos en el punto dos, se obtuvieron también dichos marcos de datos en R.

Pregunta seis

A continuación se van a crear las bases de datos con las cuales se van a hacer los ejercicios de práctica.

# Borrado de los datos

rm(list = ls())

# Base de datos uno

A <- c('a', 'b', 'c')
B <- c('t', 'u', 'v')
C <- 1:3

bd1 <- data.frame(A, B, C)
bd1

##   A B C
## 1 a t 1
## 2 b u 2
## 3 c v 3

# Base de datos dos

rm(A); rm(B); rm(C)

A <- c('b', 'c', 'd')
B <- c('u', 'v', 'w')
D <- 3:1

bd2 <- data.frame(A, B, D)
bd2

##   A B D
## 1 b u 3
## 2 c v 2
## 3 d w 1

• Primer base de datos presentada.

bd3 <- left_join(bd1, bd2, by = c('A', 'B'))
bd3

##   A B C  D
## 1 a t 1 NA
## 2 b u 2  3
## 3 c v 3  2

• Segunda base de datos

bd4 <- right_join(bd1, bd2, by = c('A', 'B'))
bd4

##   A B  C D
## 1 b u  2 3
## 2 c v  3 2
## 3 d w NA 1

• Elementos comunes

bd5 <- inner_join(bd1, bd2, by = c('A', 'B'))
bd5

##   A B C D
## 1 b u 2 3
## 2 c v 3 2

Pregunta siete

El acrónimo r4ds significa en inglés R for Data Science. El número cuatro se asocia a la palabra for dado que son palabras homófonas en este idioma. La traducción de este título en castellano es R para ciencia de datos. Los autores de este libro son Hadley Wickham y Garrett Grolemund.

Punto uno

A continuación se presenta la simulación de las observaciones para el modelo normal.

n <- 5
varianza <- 16
set.seed(202210)
xi <- runif(n = n, min = -5, max = 6)
media <- 4 - 6 * xi
set.seed(202210)
y <- rnorm(n = 5, mean = media, sd = sqrt(varianza))
datos <- data.frame(xi, y)
datos %>% kable(digits = 2, col.names = c('x', 'y'))

El diagrama de dispersión es el siguiente:

f1 <- ggplot(data = datos, aes(x = xi, y = y)) +
  geom_point() +
  ggtitle('Diagrama de dispersión para x vs. y') +
  theme_light()
f1 %>% ggplotly()

Nótese que el vector de parámetros de este modelo es:

\[ \boldsymbol{\theta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \sigma^2 \end{pmatrix} = \begin{pmatrix} 4 \\ -5 \\ 16 \end{pmatrix} \]

Ahora se va proceder a ajustar un modelo de regresión lineal de la forma:

\[ Y_i = \beta_0 + \beta_1 X_{1i} + \varepsilon_i, \ \ \ \varepsilon_i \stackrel{i.i.d}{\sim} N(0, \ \sigma^2), \ \ \ i = 1, 2, 3, 4, 5 \]

modelo <- lm(y ~ xi, data = datos)
modelo %>% stargazer(type = 'html', header = FALSE)

De acuerdo con este resultado, se tiene que el vector de parámetros estimados es:

\[ \hat{\boldsymbol{\theta}} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ \hat{\sigma^2} \end{pmatrix} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ RMSE^2 \end{pmatrix} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ RMSE^2 \end{pmatrix} = \begin{pmatrix} 0.900 \\ -7.002 \\ 2.252^2 \end{pmatrix} \begin{pmatrix} 0.900 \\ -7.002 \\ 5.07 \end{pmatrix} \]

Nótese que los valores estimados no son tan correctos respecto a los valores reales con los cuales se corrió la simulación, lo cual se debe probablemente al hecho de que solo se emplearon cinco observaciones para simular los datos.

Se va a probar de nuevo considerando ahora \(n = 300\) datos. En esta ocasión se omite la presentación de los datos simulados por la gran extensión, pero de todos modos se presenta el diagrama de dispersión.

rm(list = ls())

# Simulación de los datos

n <- 300
varianza <- 16
set.seed(202210)
xi <- runif(n = n, min = -5, max = 6)
media <- 4 - 6 * xi
set.seed(202210)
y <- rnorm(n = n, mean = media, sd = sqrt(varianza))
datos <- data.frame(xi, y)

# Gráfico de dispersión

f2 <- ggplot(data = datos, aes(x = xi, y = y)) +
  geom_point() +
  ggtitle('Diagrama de dispersión para x vs. y') +
  theme_light()
f2 %>% ggplotly()

modelo <- lm(y ~ xi, data = datos)
modelo %>% stargazer(type = 'html', header = FALSE)

En este caso se tiene que:

\[ \hat{\boldsymbol{\theta}} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ \hat{\sigma^2} \end{pmatrix} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ RMSE^2 \end{pmatrix} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\beta_1} \\ RMSE^2 \end{pmatrix} = \begin{pmatrix} 3.867 \\ -5.935 \\ 4.042^2 \end{pmatrix} \begin{pmatrix} 3.867 \\ -5.935 \\ 16.34 \end{pmatrix} \]

Se puede evidenciar que los parámetros estimados son mucho más adecuados, ya que se aproximan mucho mejor a los valores reales, evidenciando que una mayor cantidad de datos deriva en una mejor estimación.

Punto dos

El siguiente chunk muestra el bloque de código que debe ser corregido.

ni <- 50
G <- &&&
nobs <- ni * &&&                      # Numero total de observaciones
grupo <- factor(rep(x=1:G, each=&&&)) # Para crear la variable grupal
obs <- rep(x=1:ni, times=G)           # Para identificar las obs por grupo
x <- runif(n=nobs, min=&&&, max=&&&)  # La covariable
b0 <- rnorm(n=G, mean=&&&, sd=&&&)    # El Intercepto aleatorio
b0 <- rep(x=b0, each=ni)              # El intercepto aleatorio pero repetido
media <- &&& - 6 * x + &&&            # La media
&&& <- rnorm(n=nobs, mean=media, sd=&&&) # La variable respuesta
datos <- data.frame(grupo, obs, b0, x, &&&) # Organizando el dataframe

Y al hacer las correcciones correspondientes, se tiene que:

ni <- 50
G <- 10
nobs <- ni * G                           # Numero total de observaciones
grupo <- factor(rep(x=1:G, each=ni))     # Para crear la variable grupal
obs <- rep(x=1:ni, times=G)              # Para identificar las obs por grupo
set.seed(202210)
x <- runif(n=nobs, min=-5, max=6)        # La covariable
set.seed(202210)
b0 <- rnorm(n=G, mean=0, sd=sqrt(625))   # El Intercepto aleatorio
b0 <- rep(x=b0, each=ni)                 # El intercepto aleatorio pero repetido
media <- 4 - 6 * x + b0                  # La media
y <- rnorm(n=nobs, mean=media, sd= sqrt(16)) # La variable respuesta
datos <- data.frame(grupo, obs, b0, x, y) # Organizando el dataframe

f3 <- ggplot(datos, aes(x, y, color=grupo) ) + 
  geom_point() + 
  labs(colour="Grupo/Cluster")
f3 %>% ggplotly()

Puntos cuatro y cinco

Se puede observar que todas las agrupaciones tienen una tendencia lineal negativa, lo cual ocurre gracias al hecho que la pendiente no tiene un efecto aleatorio, pero sí lo tiene el intercepto, lo que se refleja en el hecho que para cada agrupación su tendencia se da en distintos niveles.

Punto seis

Los elementos del vector de parámetros del modelo con intercepto aleatorio es:

\[ \boldsymbol{\theta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \sigma^2 \\ \sigma_{b_0}^2 \end{pmatrix} = \begin{pmatrix} 4 \\ -6 \\ 16 \\ 625 \end{pmatrix} \]

Punto siete

El código a corregir es:

fit <- lmer(&&& ~ x + (1 | &&&), data=&&&)
summary(fit)

fit <- lmer(y ~ x + (1 | grupo), data = datos)
summary(fit)

## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | grupo)
##    Data: datos
## 
## REML criterion at convergence: 2840.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.01285 -0.62888  0.04132  0.65680  2.76042 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  grupo    (Intercept) 672.71   25.937  
##  Residual              14.83    3.851  
## Number of obs: 500, groups:  grupo, 10
## 
## Fixed effects:
##             Estimate Std. Error  t value
## (Intercept)  0.76828    8.20376    0.094
## x           -6.07390    0.05505 -110.334
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.003

Estilizado, luce como sigue:

En este caso,

\[ \hat{\boldsymbol{\theta}} = \begin{pmatrix} \hat{\beta}_0 \\ \hat{\beta}_1 \\ \hat{\sigma}^2 \\ \hat{\sigma}_{b_0}^2 \end{pmatrix} = \begin{pmatrix} 0.76828 \\ -6.07390 \\ 14.83 \\ 672.71 \end{pmatrix} \]

Se puede evidenciar que todas las estimaciones son cercanas a los valores reales empleados en la simulación, excecto el del intercepto \(\hat{\beta}_0\)

Una medida de cercanía entre valor es el error absoluto relativo, el cual se calcula como sigue:

\[ EAR = \frac{| \text{Valor verdado - valor estimado} |}{\text{valor verdadero}} \times 100 \% \] Y se pueden calcular rápidamente con ayuda de R como se muestra enseguida:

reales <- c(4, -6, 16, 625)
estimados <- c(0.76828, -6.07390, 14.83, 672.71)
ear <- abs(reales - estimados) * 100 / reales
nombres <- c('Intercepto fijo', 'Pendiente fija',
             'Varianza de la respuesta', 'Varianza del intercepto aleatorio')
errores <- data.frame(reales, estimados, ear)
row.names(errores) <- nombres
errores %>% kable(digits = 2,
                  col.names = c('Valor real', 'Valor estimado',
                                'EAR [%]'))

Pregunta diez

cbind(Verdadero = unique(b0), Estimado = ranef(fit)$grupo) %>% kable(digits = 3)

Luce como si los resultados fueran buenos.

Punto uno

El nombre de la base de datos que está siendo usada en el video es hsb y se encuentra en el paquete merTools.

Punto dos

Al emplear la función str() con esta base de datos, se puede observar que esta contiene ocho variables, a saber:

• schid. Es una variable tipo carácter que contiene un identificador para cada escuela considerada en esta base de datos.
• minority. Es una variable dicotómica que indica si el estudiante pertenece o no a una minoría. Si toma el valor de cero es porque el estudiante no hace parte de ninguna minoría (el estudiante es blanco) y toma el valor de uno en el caso contrario.
• female. Es una variable dicotómica asociada con el género del estudiante, la cual toma el valor de cero si el estudiante es de género masculino y toma el valor de uno cuando la estudiante es de género femenino.
• ses. Es una variable numérica que mide el estatus socioeconómico del estudiante. . mathach. Es una variable numérica que indica el logro en un test de matemáticas por parte del estudiante.
• size. Es una variable numérica común a cada una de las escuelas.
• schtype. Es una variable dicotómica que indica el tipo de financiación de la escuela, por lo que toma el valor de cero cuando es una institución pública y uno si es privada.
• meanses. Es un índice numérico continuo que promedia la situación socioecónomica de los estudiantes de cada escuela, de tal suerte que esta característica va a ser común para cada escuela.

Punto tres

La variable respuesta (\(y\)) que se va a tratar de predecir usando a las demás como covariables es el logro en matemáticas, es decir, la característica denominada como mathach.

Punto cuatro

El modelo uno es un modelo de intercepto aleatorio de acuerdo con la variable schid, que es la variable que está siendo empleanda como variable de agrupación.

Punto cinco

El parámetro lógico REML que se emplea en la función lme4::lmer() sirve para emplear el método de estimación que se va emplear para los parámetros de interés de acuerdo con el modelo a ajustar. Así, si este parámetro es FALSE, se va a calcular \(\hat{\boldsymbol{\theta}}\) usando máxima verosimilitud, mientras que si este es TRUE, se utiliza máxima verosimilitud parcial o restringida (el nombre viene de restricted o reduced maximum likelihood).

Punto seis

El vector de parámetros \(\boldsymbol{\theta}\) asociado al modelo uno solo va a tener tres parámetros, ya que están ajustando modelos que solo tienen intercepto aleatorio de acuerdo con la escuela. En otras palabras, el vector de parámetros \(\boldsymbol{\theta}\) es:

\[ \boldsymbol{\theta} = \begin{pmatrix} \beta_0 \\ \sigma^2 \\ \sigma_{b_0}^2 \end{pmatrix} \]

Punto siete

El coeficiente de correlación interclase es una métrica que se emplea para medir la relación que existe entre dos variables de distintas clases. En este caso, por ejemplo, podría interesar la correlación que existe en los puntajes de logro en matemáticas entre estudiantes de una escuela A y los estudiantes de una escuela B. Esta métrica resulta, pues, útil cuando se están considerando datos que están agrupados o clusterizados, como es el caso en esta situación según la característica schid. Así, sea \(\sigma_{b_0}^2\) la varianza asociada a los grupos y sea \(\sigma_y^2\) la varianza general, entonces esta métrica se puede calcular empleando la siguiente ecución.

\[ ICC = \frac{\sigma_{b_0}^2}{\sigma_{b_0}^2 + \sigma_y^2} \]

Punto ocho

Usando el modelo uno, se tiene que:

\[ \hat{\boldsymbol{\theta}} = \begin{pmatrix} \hat{\beta_0} \\ \hat{\sigma}^2 \\ \hat{\sigma}_{b_0}^2 \end{pmatrix} = \begin{pmatrix} 12.6371 \\ 39.148 \\ 8.553 \end{pmatrix} \]

Con esto presente, se tiene que de acuerdo con este modelo, el coeficiente de correlación interclase es:

\[ ICC = \frac{\sigma_{b_0}^2}{\sigma_{b_0}^2 + \sigma_y^2} \approx \frac{8.553}{8.553 + 39.148} \approx 0.18 \]

Punto uno

Considerar el siguiente modelo normal mixto con pendiente aleatoria \(b_{1i}\):

\[ \begin{cases} y_i \sim N(\mu_{ij}, \ \sigma_y^2) \\ \mu_{ij} = 4 - 5x_{ij} + b_{1i} x_{ij} \\ \sigma_y^2 = 4 \\ b_1 \sim N(0, \ \sigma_{b1}^2 = 1000) \\ x_{ij} \sim Beta (\alpha = 2, \ \beta = 4) \end{cases} \]

Escribir los cuatro elementos del vector de parámetros \(\boldsymbol{\Theta}\) del modelo.

En este caso se tiene que:

\[ \boldsymbol{\Theta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \sigma_y^2 \\ \sigma_{b1}^2 \end{pmatrix} = \begin{pmatrix} 4 \\ -5 \\ 4 \\ 1000 \end{pmatrix} \]

Punto dos

Completar el siguiente código de \(\color{blue}{\texttt{R}}\) para simular \(n_i = 50\) observaciones para \(G = 10\) grupos, es decir, en total \(500\) observaciones. Cambiar los símbolos &&& por el código correcto para simular los datos solicitados.

El código dado es el siguiente:

ni <- &&&
G <- &&&
nobs <- &&& * &&&                          # Numero total de observaciones
grupo <- factor(rep(x=1:G, &&&=ni))        # Para crear la variable grupal
obs <- rep(x=1:ni, times=G)                # Para identificar las obs por grupo
&&& <- rbeta(n=nobs, &&&=2, &&&=&&&)       # La covariable x
b1 <- rnorm(n=&&&, mean=&&&, sd=sqrt(&&&)) # La pendiente aleatoria
&&& <- rep(x=b1, &&&=&&&)                  # La pendiente repetida
media <- &&& + &&& * x + b1 * x            # La media
y <- &&&(n=nobs, mean=&&&, sd=sqrt(&&&))   # La variable respuesta
datos <- &&&(grupo, obs, b1, x, &&&)       # Organizando el dataframe

Ahora, se van a hacer dos correcciones: se van a cambiar los signos &&& por el código correspondiente, y se va a crear una función llamada datos_simulados() que facilite la generación de los datos, teniendo como parámetro a semilla, la cual permite generar siempre los mismos datos. Esta función va a regresar un marco de datos con los datos simulados.

datos_simulados <- function(semilla = semilla, ni = 50, G = 10, sigma_2_b1 = 1000 ) {
  ni <- ni
  G <- G
  nobs <- ni * G
  grupo <- factor(rep(x = 1:G, each = ni))
  obs <- rep(x = 1:ni, times = G)
  set.seed(semilla)
  x <- rbeta(n = obs, shape1 = 2, shape2 = 4)
  set.seed(semilla)
  b1 <- rnorm(n = G, mean = 0, sd = sqrt(sigma_2_b1))
  b1 <- rep(x = b1, each = ni)
  media <- 4 - 5 * x + b1 * x
  y <- rnorm(n = obs, mean = media, sd = sqrt(4))
  datos <- data.frame(grupo, obs, b1, x, y)
  return(datos)
}

Es importante notar que la función anterior no es generalizable, es decir, solo sirve para simular datos bajo la estructura indicada en el punto uno, bajo una semilla de generación de números pseudoaleatorios definida por el usuario para asegurar la replicabilidad de los datos, por lo que si se modifica algún parámetro, esta debe cambiar (salvo para \(\sigma_{b1}^2\) de tal forma que se pueda reciclar esta función en el punto seis). Lo único en lo que esta función es flexible es en la cantidad de datos a simular, si bien sus valores predeterminados son \(n_i = 50\) y \(G = 10\).

Ahora se procede a chequear los resultados de esta función usando a \(1037670103\) como semilla de generación. Se mostrarán las primeras cinco y las últimas cinco observaciones simuladas.

datos <- datos_simulados(semilla = 20221010)
datos[c(1:5, 496:500), ] |> knitr::kable(digits = 2)

Punto tres

Ilustrar los datos en un diagrama de dispersión, coloreando los datos según la agrupación a la que pertenecen.

ggplot(datos, aes(x, y, color = grupo)) + 
  geom_point(size = 1, pch = 20) + 
  labs(colour = "Agrupación")

Como se puede observar, todas las agrupaciones tienen un intercepto común, que se aproxima a cero (como se espera a partir de la simulación). No obstante, se evidencia que cada agrupación tiene una pendiente diferente, la cual es determinada por el efecto aleatorio \(b_1\).

Punto cuatro

Encontrar las estimaciones \(\hat{\boldsymbol{\Theta}}\) del vector \(\boldsymbol{\Theta}\). ¿Qué tan cerca está la estimación del vector real?

Se va a usar la función lmer() del paquete lme4 para poder obtener el vector de estimaciones \(\hat{\boldsymbol{\Theta}}\):

modelo <- lmer(y ~ x + (x - 1 | grupo), REML = FALSE, data = datos)
modelo |> summary()

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: y ~ x + (x - 1 | grupo)
##    Data: datos
## 
##      AIC      BIC   logLik deviance df.resid 
##   2208.0   2224.9  -1100.0   2200.0      496 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8415 -0.6677 -0.0125  0.5854  3.4150 
## 
## Random effects:
##  Groups   Name Variance Std.Dev.
##  grupo    x    408.004  20.199  
##  Residual        4.188   2.046  
## Number of obs: 500, groups:  grupo, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    3.929      0.198  19.839
## x             -3.133      6.410  -0.489
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.074

De acuerdo con la salida de \(\color{blue}{\texttt{R}}\), se tiene que:

\[ \hat{\boldsymbol{\Theta}} = \begin{pmatrix} \hat{\beta}_0 \\ \hat{\beta}_1 \\ \hat{\sigma}_y^2 \\ \hat{\sigma}_{b1}^2 \end{pmatrix} = \begin{pmatrix} 3.923 \\ -3.133 \\ 4.188 \\ 408.004 \end{pmatrix} \]

Se puede observar que las estimaciones asociadas a los efectos fijos son adecuadas, pues no difieren significativamente de los valores reales. Algo similar sucede con la estimación de la varianza del efecto aleatorio asociado a la respuesta, pero no pasa igual con la varianza estimada para el efecto aleatorio de la pendiente, que es menos de la mitad de la real.

Punto cinco

Encontrar las predicciones de los \(b_1\) y compararlas con los valores reales.

Se va a calcular el error absoluto relativo así:

\[ \text{Error absoluto relativo} = \frac{|\text{Valor real} - \text{Valor estimado}|}{\text{Valor real}} \]

estimaciones <- ranef(modelo)$grupo
reales <- datos$b1 |> unique()
errores <- abs(estimaciones - reales)/reales
grupo <- 1:10
comparaciones <- data.frame(grupo, estimaciones, reales, errores)
comparaciones |> knitr::kable(digits = 2,
                              col.names = c('Grupo', 'Valor real', 
                                            'Valor estimado',
                                            'Error absoluto relativo'))

Mientras más cercano sea el error absoluto relativo a cero, mejor. Y como se observa, para la mayoría de las variables esto sucede con errores relativos del orden del 10 % o menos, excepto para el valor estimado asociado al grupo seis, con una sobrestimación del 20 %, pero es peor aún para el grupo tres, ya que en este caso se subestimó el valor hasta su tercera parte aproximadamente.

Punto seis

Simular de nuevo el marco de datos, llamado ahora datos2 pero con \(\sigma_{b1}^2 = 0\) y dibujar los datos. ¿Qué se observa ahora con las nubes de puntos?

# Simulación de los datos
datos2 <- datos_simulados(semilla = 20221010, sigma_2_b1 = 0)

# Graficación de los datos
ggplot(datos2, aes(x, y, color = grupo)) + 
  geom_point(size = 1, pch = 20) + 
  labs(colour = "Agrupación")

En este caso se observa una nube de puntos con una apariencia más aleatoria, por así decirlo, que la anterior, ya que en este caso, como la varianza es nula, entonces la variable aleatoria en realidad es degenerada puesto que es una constante, de tal forma que todos los grupos tendrán la misma pendiente de aproximadamente cinco, siendo gobernados ahora por la variable \(x\). Lo anterior implica que la pendiente es \(-5\), y esto se refleja en la tendencia negativa que tienen los puntos.

Punto uno

Se puede notar que se tiene un modelo jerárquico de pendiente e intercepto aleatorios. Así, se van a tener siete elementos en el vector de parámetros, a saber:

\[ \boldsymbol{\Theta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \sigma_y^2 \\ \sigma_{b0}^2 \\ \sigma_{b1}^2 \\ \sigma_{b01} \end{pmatrix} = \begin{pmatrix} 4 \\ -5 \\ 4 \\ 40 \\ 50 \\ 3 \end{pmatrix} \]

Punto dos

A continuación se muestra el bloque de código con los elementos que se deben reemplazar.

ni <- &&&
G <- &&&
nobs <- &&& * &&& # Numero total de observaciones
grupo <- factor(rep(x=1:G, &&&=ni)) # Para crear la variable grupal
obs <- rep(x=1:ni, &&&=G) # Para identificar las obs por grupo
x <- &&& # La covariable
library(MASS) # Libreria para simular obs de Normal bivariada
Sigma <- matrix(c(40, 3, # Matriz de var-cov
                  3, 50), ncol=2, nrow=2)
b <- mvrnorm(n=G, mu=rep(0, 2), Sigma=Sigma) # Simulando b0 y b1
b <- apply(b, MARGIN=2, function(c) rep(c, each=ni)) # Replicando
b0 <- as.vector(b[, 1]) # Separando los b0
b1 <- as.vector(b[, 2]) # Separando los b1
media <- &&& + &&& * x + b0 + b1 * x # La media
y <- rnorm(n=&&&, mean=&&&, sd=&&&) # La variable respuesta
datos <- &&&(grupo, obs, &&&, &&&, &&&, &&&) # El dataframe

De manera que al hacer las correcciones pertinentes, se obtienen los siguientes resultados:

# Parámetros iniciales generales
ni <- 50                                     # Observaciones por grupo
G <- 10                                      # Número de grupos
nobs <- ni * G                               # Cantidad total de observaciones
grupo <- rep(x = 1:G, each = ni) %>% 
  factor()                                   # Variable de agrupación
obs <- rep(x = 1:ni, times = G)              # Identificador del número de la ~
                                             # observación dentro de cada grupo

# Simulación de la covariable
set.seed(202210)                             # Semilla para garantizar repeti~
                                             # bilidad
xi <- runif(n = nobs, min = 0, max = 5)      # Xi ~ Uniforme(0, 5)

# Matriz de varianzas y covarianzas
Sigma <- matrix(c(40, 3,
                  3, 50),
                ncol = 2, nrow = 2)

# Simulación de la pendiente y el intercepto aleatorios
set.seed(202210)                             # Semilla para garantizar repeti~
                                             # bilidad
b <- mvrnorm(n = G, mu = c(0, 0),
             Sigma = Sigma)                  # Simulación de b0 y b1
b <- apply(b, MARGIN = 2,
           function(c)
             rep(c, each = ni))              # Replicando ni veces en c/ grupo
b0 <- as.vector(b[, 1])                      # Separando los b0
b1 <- as.vector(b[, 2])                      # Separando los b1

# Variable de respuesta
media <- (4 + b0) + (b1 - 5) * xi            # Ecuación de la media
sd_y <- 4                                    # Varianza de la respuesta
set.seed(202210)                             # Semilla para garantizar repeti~
                                             # bilidad
y <- rnorm(n = nobs, mean = media, sd = sd_y)# yij ~ N(media, 4)

# Marco de datos
datos <- data.frame(grupo, obs, b0, b1, xi, y)

A continuación se expone cómo lucen las primeras y las últimas cinco observaciones de este marco de datos:

exposicion <- rbind(head(datos, 5), tail(datos, 5))
nombres <- c('Grupo', 'Observación en el grupo', 
             'Intercepto aleatorio', 'Pendiente aleatoria',
             'X', 'Y')
exposicion %>% kable(digits = 2, col.names = nombres)

Punto tres

f4 <- ggplot(datos, aes(xi, y, color = grupo)) + 
  geom_point(size = 1, pch = 20) + 
  labs(colour = "Agrupación") +
  ggtitle('Datos simulados',
          subtitle = 'Para un modelo mixto de pendiente e intercepto aleatorios') +
  xlab('X') +
  ylab('Y') +
  scale_color_npg() +
  theme_light()
f4 %>% ggplotly()

Punto cuatro

modelo_c <- lmer(y ~ xi + (xi + 1 | grupo),
                 REML = FALSE, data = datos)
modelo_c %>% summary()

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: y ~ xi + (xi + 1 | grupo)
##    Data: datos
## 
##      AIC      BIC   logLik deviance df.resid 
##   2874.0   2899.3  -1431.0   2862.0      494 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8583 -0.6510  0.0240  0.6691  2.9576 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  grupo    (Intercept) 25.96    5.095        
##           xi          49.65    7.046    0.07
##  Residual             14.57    3.816        
## Number of obs: 500, groups:  grupo, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    4.418      1.648   2.681
## xi            -6.048      2.231  -2.711
## 
## Correlation of Fixed Effects:
##    (Intr)
## xi 0.058

\[ \hat{\boldsymbol{\Theta}} = \begin{pmatrix} \hat{\beta}_0 \\ \hat{\beta}_1 \\ \hat{\sigma}_y^2 \\ \hat{\sigma}_{b0}^2 \\ \hat{\sigma}_{b1}^2 \\ \hat{\sigma}_{b01} \end{pmatrix} = \begin{pmatrix} 4.418 \\ -6.048 \\ 14.57 \\ 25.96 \\ 49.65 \\ 2.513 \end{pmatrix} \]

• Pregunta. ¿Es acertado el cálculo que hice para la covarianza? Correlación * sd_b0 * sd_b1

Se puede evidenciar que \(\beta_0, \ \beta_1, \ \sigma_{b1}^2\) y \(\sigma_{b01}\) tienes estimaciones cercanas a sus valores reales de acuerdo con la simulación realizada. Este no es el caso para \(\sigma_y^2\) y para \(\sigma_{b0}^2\)

Pregunta cinco

# Interceptos reales
b0_real <- datos$b0 %>% unique()
b0_real <- b0_real + 4

# Pendientes reales
b1_real <- datos$b1 %>% unique()
b1_real <- b1_real - 5

# Estimados
estimados <- (modelo_c %>% coef())$grupo
estimados %<>% janitor::clean_names()
int_est <- estimados$intercept
pen_est <- estimados$xi

# Errores realtivos
ear_b0 <- abs((b0_real - int_est) / b0_real) * 100
ear_b1 <- abs((b1_real - pen_est) / b1_real) * 100

# Comparaciones
res <- data.frame(b0_real, int_est, ear_b0,
             pen_est, estimados$xi, ear_b1)

# Exposición
nombres2 <- c('Intercepto real', 'Intercepto estimado', 'EAR del intercepto [%]',
              'Pendiente real', 'Pendiente estimada', 'EAR de la pendiente [%]')
res %>% kable(digits = 2, col.names = nombres2)

Punto uno

Se puede observar en el panel dos que las rectas de regresión ajustadas a cada agrupación son paralelas, lo cual indica de inmediato que se trata de un modelo de intercepto aleatorio. A continuación se desarrolla su ajuste:

modelo1 <- lmer(Reaction ~ Days + (1 | Subject),
                REML = TRUE, data = sleepstudy)
modelo1 %>% summary()

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

Y se replica la gráfica enseguida:

coefs1 <- ((modelo1 %>% coef)$Subject %>% janitor::clean_names())$intercept
fit <- predict(modelo1)

f5 <- ggplot(sleepstudy, aes(x = Days, y = Reaction, color = Subject)) + 
        geom_point(size = 1.5, pch = 20) + 
        labs(colour = "Agrupación") +
        geom_line(y = fit) +
        ggtitle('Reación según el número de días transcurridos por sujeto',
                subtitle = 'Modelo de intercepto aleatorio') +
        xlab('Días') +
        ylab('Tiempo de reacción [min]') +
        theme_light()
f5 %>% plotly::ggplotly()

A la luz de los resultados, este modelo no parece ser muy acertado ya que si bien parece razonable incluir el intercepto aleatorio, no todos los individuos tienen la misma tasa de cambio en cuanto al tiempo de reacción a medida que pasan los días.

Punto dos

Se puede observar en el panel tres que las rectas de regresión ajustadas a cada agrupación parten del mismo punto, lo cual indica de inmediato que se trata de un modelo de pendiente aleatoria. A continuación se desarrolla su ajuste:

modelo2 <- lmer(Reaction ~ Days + (Days - 1 | Subject),
                REML = TRUE, data = sleepstudy)
modelo2 %>% summary()

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days - 1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1766.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5104 -0.5588  0.0541  0.6244  4.6022 
## 
## Random effects:
##  Groups   Name Variance Std.Dev.
##  Subject  Days  52.71    7.26   
##  Residual      842.03   29.02   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   251.41       4.02  62.539
## Days           10.47       1.87   5.599
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.340

Y se replica la gráfica enseguida:

fit2 <- predict(modelo2)

f6 <- ggplot(sleepstudy, aes(x = Days, y = Reaction, color = Subject)) + 
        geom_point(size = 1.5, pch = 20) + 
        labs(colour = "Agrupación") +
        geom_line(y = fit2) +
        ggtitle('Reación según el número de días transcurridos por sujeto',
                subtitle = 'Modelo de pendiente aleatoria') +
        xlab('Días') +
        ylab('Tiempo de reacción [min]') +
        theme_light()
f6 %>% plotly::ggplotly()

Se puede pensar que este modelo es aún más inadecuado que el anterior, ya que si bien la tasa de cambio no es la misma para todos los sujetos, sí son aproximadamente parecidas entre sí, mientras que de acuerdo con la gráfica anterior no resulta claro que los interceptos sean los mismos y las tasas de cambio del tiempo de reacción tan dispares entre sí.