Este ejercicio trata sobre cuando el supuesto de varianzas iguales no se cumple en el modelo de regresión lineal, una de las alternativas es una regresión ponderada, supongamos que tenemos los siguientes datos:
cant <- c(5, 6, 7, 8, 9)
rend <- c(2, 2.8, 3, 1, 0.5)
n <- c(3, 7, 2, 8, 4)
data.frame(Cantidad = cant, Rendimiento = rend, Replicas = n)
## Cantidad Rendimiento Replicas
## 1 5 2.0 3
## 2 6 2.8 7
## 3 7 3.0 2
## 4 8 1.0 8
## 5 9 0.5 4
Donde la variable explicativa corresponde a la cantidad (realmente no se de que), notemos que los datos se encuentran agregados en una tabla de contingencia, con el rendimiento promedio como variable explicada. Siguiendo algunas de las instrucciónes del modelo, ajustaremos un modelo polinomial cuadrático, lo cual también nos ayudara a mostrar gráficamente algunos puntos. De tal manera que el modelo a ajustar sería:
\[ \bar{y}=\beta_0+\beta_1 x+\beta_2 x^2 \]
Sin embargo como el número de réplicas es diferente en cada grupo, y no tenemos los datos desagregados, el supuesto de homocedasticidad no se cumple. Así que para corregir esto ajustamos el modelo ponderado:
\[ \sqrt{n_i}\bar{y}=\sqrt{n_i}\beta_0+\sqrt{n_i}\beta_1 x+\sqrt{n_i}\beta_2 x^2 \]
Ajustamos ambos modelos:
mod_ols <- lm(rend ~ cant + I(cant^2)) #Modelo de mínimos cuadrados ordinarios
mod_pond <- lm(I(sqrt(n) * rend) ~ 0 + sqrt(n) + I(sqrt(n) * cant) + I(sqrt(n) *
cant^2)) #Modelo ponderado
Notemos las diferencias entre los coeficientes, es decir que son modelos diferentes:
coef(mod_ols)
## (Intercept) cant I(cant^2)
## -10.8943 4.3200 -0.3429
coef(mod_pond)
## sqrt(n) I(sqrt(n) * cant) I(sqrt(n) * cant^2)
## -6.4713 3.0562 -0.2583
Ahora bien, dentro de R tenemos otra forma de de incluir los pesos en el modelo (weights), esto es directamente en la función de R lm, lo haremos de dos maneras, la primera utilizando la probabilidad observada en cada grupo y la segunda utilizando el número de réplicas en cada grupo, notemos que para ambos casos los coeficientes son los mismos que el modelo ponderado:
w <- n/sum(n)
mod_w <- lm(rend ~ cant + I(cant^2), weights = w) #Probabilidades
mod_n <- lm(rend ~ cant + I(cant^2), weights = n) #Réplicas
coef(mod_w)
## (Intercept) cant I(cant^2)
## -6.4713 3.0562 -0.2583
coef(mod_n)
## (Intercept) cant I(cant^2)
## -6.4713 3.0562 -0.2583
Gracias a que es un modelo polinomial con una sola variable regresora podemos observar la diferencia entre los modelos gráficamente:
x_seq <- data.frame(cant = seq(from = 4, to = 10, by = 0.05))
plot(data.frame(cantidad = x_seq, rendimiento = predict(mod_ols, x_seq)), type = "l",
main = expression(bar(y) == beta[0] + beta[1] * x + beta[2] * x^2), col = 2,
xlab = "Cantidad")
abline(h = 0, lty = 2)
lines(data.frame(cantidad = x_seq, rendimiento = predict(mod_w, x_seq)), col = 4)
legend("bottomleft", bty = "n", legend = c("OLS", "Ponderado"), col = c(2, 4),
lty = 1)
Resulta evidente que si bien las curvas son similares, es claro que no son iguales. Cabe notar que el modelo no es bueno, para empezar por que con una cantidad de 10, el rendimiento se vuelve negativo (aunque no sabemos si el rendimiento puede ser negativo, pero parece poco probable)
Ahora pasemos a la siguiente cuestión, tenemos tres modelos con coeficientes esetimados iguales, pero ¿En realidad son los mismos modelos? Observemos más detalladamente los resultados:
summary(mod_pond)
##
## Call:
## lm(formula = I(sqrt(n) * rend) ~ 0 + sqrt(n) + I(sqrt(n) * cant) +
## I(sqrt(n) * cant^2))
##
## Residuals:
## 1 2 3 4 5
## -0.611 0.615 1.038 -1.267 0.774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## sqrt(n) -6.471 9.763 -0.66 0.58
## I(sqrt(n) * cant) 3.056 2.846 1.07 0.40
## I(sqrt(n) * cant^2) -0.258 0.202 -1.28 0.33
##
## Residual standard error: 1.42 on 2 degrees of freedom
## Multiple R-squared: 0.957, Adjusted R-squared: 0.893
## F-statistic: 14.9 on 3 and 2 DF, p-value: 0.0637
summary(mod_w)
##
## Call:
## lm(formula = rend ~ cant + I(cant^2), weights = w)
##
## Weighted Residuals:
## 1 2 3 4 5
## -0.125 0.125 0.212 -0.259 0.158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.471 9.763 -0.66 0.58
## cant 3.056 2.846 1.07 0.40
## I(cant^2) -0.258 0.202 -1.28 0.33
##
## Residual standard error: 0.29 on 2 degrees of freedom
## Multiple R-squared: 0.815, Adjusted R-squared: 0.63
## F-statistic: 4.4 on 2 and 2 DF, p-value: 0.185
summary(mod_n)
##
## Call:
## lm(formula = rend ~ cant + I(cant^2), weights = n)
##
## Weighted Residuals:
## 1 2 3 4 5
## -0.611 0.615 1.038 -1.267 0.774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.471 9.763 -0.66 0.58
## cant 3.056 2.846 1.07 0.40
## I(cant^2) -0.258 0.202 -1.28 0.33
##
## Residual standard error: 1.42 on 2 degrees of freedom
## Multiple R-squared: 0.815, Adjusted R-squared: 0.63
## F-statistic: 4.4 on 2 and 2 DF, p-value: 0.185
Resulta que tanto las estimaciones de los coeficientes como las desviaciones de estos entre los modelos son iguales. Sin embargo entre el modelo ponderado y los modelos que utilizaron la función weights difiere el estadístico F, sobre la significacia del modelo, aunque en ambos casos la concusión es la misma: el modelo no resulta significativo. Pero entre los modelos que utilizaron la función weights los residuales difieren; cuando se utilizó el vector de réplicas los residuales son idénticos al modelo ponderado visto en clase. Veamos lo que dice la ayuda de R al respecto:
weights: an optional vector of weights to be used in the fitting process. Should be NULL or a numeric vector. If non-NULL, weighted least squares is used with weights weights (that is, minimizing sum(w*e^2)); otherwise ordinary least squares is used. See also ‘Details’...
Non-NULL weights can be used to indicate that different observations have different variances (with the values in weights being inversely proportional to the variances); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations (including the case that there are w_i observations equal to y_i and the data have been summarized).
En este caso resulta que lo que debemos de utilizar es en realidad, dado que desconocemos la varianza dentro de cada grupo, debemos de utilizar el vector de réplicas, nos olvidaremos entonces del otro modelo.Tomando esto en cuenta podemos calcular los intervalos de predicción para el modelo (que dado que los residuales coinciden con el modelo visto en clase, estos intervalos coinciden también).
Esto nos lleva a la siguiente cuestión, sobre la información que se tiene, o más específicamente la información que no se tiene, es decir que en al tener la información en forma agregada perdemos información importante, para el caso de la regresión específicamente la varianza y en este caso, la varianza dentro de cada grupo. Para este caso vamos a introducir información, supuesta que nos lleve a las mismas medias pero con varianzas diferentes:
x_2 <- c(rep(cant, n))
y_2 <- c(1, 2, 3, 5.2, 4, 4.4, 2, 2, 1, 1, 4, 2, 2, 1, 0.5, 0.5, 1, 1, 2, 0,
1, 0, 0.5, 0.5)
aggregate(y_2 ~ x_2, FUN = mean)
## x_2 y_2
## 1 5 2.0
## 2 6 2.8
## 3 7 3.0
## 4 8 1.0
## 5 9 0.5
Utilizaremos este modelo y lo compararemos con el del vector de medias, notese que los coeficientes coinciden, pero veamos la diferencia tomando en cuenta los intervalos de predicción:
mod_ext2 <- lm(y_2 ~ x_2 + I(x_2^2))
coef(mod_ext2)
## (Intercept) x_2 I(x_2^2)
## -6.4713 3.0562 -0.2583
x_seq2 <- x_seq
names(x_seq2) <- "x_2"
plot(data.frame(cantidad = x_seq, Rendimiento = predict(mod_w, x_seq)), type = "l",
main = expression(bar(y) == beta[0] + beta[1] * x + beta[2] * x^2), col = 3,
ylim = c(-10, 9), lwd = 2)
abline(h = 0, lty = 2, lwd = 0.5)
lines(data.frame(cantidad = x_seq, rendimiento = predict(mod_n, x_seq, interval = "p")[,
2]), lty = 2, col = 2)
## Warning: Assuming constant prediction variance even though model fit is weighted
lines(data.frame(cantidad = x_seq, rendimiento = predict(mod_n, x_seq, interval = "p")[,
3]), lty = 2, col = 2)
## Warning: Assuming constant prediction variance even though model fit is weighted
lines(data.frame(cantidad = x_seq, rendimiento = predict(mod_ext2, x_seq2, interval = "p")[,
2]), lty = 2, col = 4)
lines(data.frame(cantidad = x_seq, rendimiento = predict(mod_ext2, x_seq2, interval = "p")[,
3]), lty = 2, col = 4)
legend("bottom", bty = "n", lty = 2, col = c(2, 4), legend = c("Réplicas", "Datos completos"))
Notemos como los intervalos del modelo con los datos completos tiene intervalos de predicción significativamente menores; ahora siguiendo la idea sobre la importancia de la varianza en el modelo, el modelo ponderado con las varianzas dentro de cada grupo debería de resultar en el mismo modelo de los datos completos, veamos si esto es así:
s_n <- aggregate(y_2 ~ x_2, FUN = var)[, 2]
s2 <- (s_n^-1) + (s_n * (n - 1))^-1
mod_var <- lm(rend ~ cant + I(cant^2), weights = I(n/s_n))
Según yo debería de ser así sin embargo probe muchas formas para introducir el vector de pesos como proporcinal inverso a la varianza, pero no funcionó, se probó con la muestral,poblacional, su inverso, considerando la muestra para cada grupo, etc.,etc., hasta que me cansé, la mejor explicación la encontré aquí, pero no funcionó.
Según yo la idea es correcta, si tenemos la varianza, y los datos agrupados,esta información es suficiente para un modelo lineal. Sin embargo, no lo pude aplicar con la función, y la verdad para ponerse a hacer algebra, no es el momento. No hay como tener los datos desagregados.