6.3 Inferencia sobre los parámetros del Modelo
Una vez evaluado el correcto ajuste del modelo utilizando las metodologías vistas anteriormente y corroborando las propiedades distribucionales de los errores y por ende, de la variable respuesta \(y\), el paso siguiente es verificar si los parámetros estimados son significativos y por ende, las covariables utilizadas para ajustar el modelo aportan significativamente para explicar o predecir a la variable de estudio \(y\).
Dada las propiedades distribucionales de los \(\beta's\), un estadístico de prueba natural para evaluar la significancia de dicho parámetro se basa en la distribución “t-student” y se describe a continuación:
\[ t=\frac{\hat{\beta}_{k}-\beta_{k}}{se\left(\hat{\beta}_{k}\right)}\sim t_{n-p} \]
Donde \(p\) es el número de parámetros del modelo y \(n\) el tamaño de la muestra de la encuesta. En este sentido, el estadístico de prueba anterior evalúa las hipótesis \(H_{0}:\beta_{k}=0\) versus la alternativa \(H_{1}:\beta_{k}\neq0\).
De las propiedades distribucionales de los \(\beta\), se puede construir un intervalo de confianza al \((1-\alpha)\times100\%\) para \(\beta_{k}\) está dado por:
\[ \hat{\beta}_{k}\pm t_{1-\frac{\alpha}{2},\,df}\,se\left(\hat{\beta}_{k}\right) \]
Donde, los grados de libertad para el intervalo (\(df\)) en una encuesta de hogares (muestras complejas) está dado por el número de conglomerados finales de la primera etapa menos el número de estratos de la etapa primaria \(\left(df=\sum_{h}a_{h}-H\right)\).
Para la aplicación de las temáticas vistas, es decir, realizar la prueba de hipótesis y los intervalos de confianza para los parámetros utilizaremos el modelo que se ha venido trabajando y aplicaremos las funciones summary.svyglm para las pruebas t y confint.svyglm para los intervalos de confianza como sigue:
survey:::summary.svyglm(mod_svy)##
## Call:
## svyglm(formula = Income ~ Expenditure + Zone + Sex + Age2, design = diseno_qwgt)
##
## Survey design:
## Called via srvyr
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.18419 62.98928 0.99 0.33
## Expenditure 1.22548 0.19798 6.19 9.5e-09 ***
## ZoneUrban 63.46000 40.09025 1.58 0.12
## SexMale 21.73256 15.78081 1.38 0.17
## Age2 0.00852 0.00557 1.53 0.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 118297)
##
## Number of Fisher Scoring iterations: 2survey:::confint.svyglm(mod_svy)| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | -62.5855 | 186.9538 |
| Expenditure | 0.8333 | 1.6176 |
| ZoneUrban | -15.9511 | 142.8711 |
| SexMale | -9.5262 | 52.9913 |
| Age2 | -0.0025 | 0.0196 |
De lo anterior se puede observar que, con una confianza del 95% el único parámetro significativo del modelo es Expenditure y ese mismo resultado lo reflejan los intervalos de confianza.
Estimación de una observación
Los modelos de regresión lineales, según (Neter et al., 1996)., son utilizado esencialmente con 2 fines, el primero es tratar de explicar la variable respuesta en términos de covariables que pueden encontrarse en la encuesta o en registros administrativos, censos, etc. Adicionalmente, también son usados para predecir valores de la variable en estudio ya sea dentro del intervalo de valores recogidos en la muestra o por fuera de dicho intervalo. Lo primero se ha abordado a lo largo de todo el capítulo y lo segundo se obtiene de la siguiente manera:
\[ \hat{E}(y_{i}\mid\boldsymbol{x}_{obs,i})=\boldsymbol{x}_{obs,i}\hat{\boldsymbol{\beta}} \]
De manera explícita, si se ajusta un modelo con 4 covariables la expresión sería:
\[ \hat{E}(y_{i}\mid\boldsymbol{x}_{obs,i})=\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\hat{\beta}_{2}x_{2i}+\hat{\beta}_{3}x_{3i}+\hat{\beta}_{4}x_{4i} \]
La varianza de la estimación se calcula de la siguiente manera:
\[ var\left(\hat{E}\left(y_{i}\mid x_{obs,i}\right)\right) = x'_{obs,i}cov\left(\hat{\beta}\right)x{}_{obs,i} \]
A continuación, se presenta cómo se realiza la estimación del valor esperado, primero se estiman los parámetros del modelo:
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 62.1842 | 62.9893 | 0.9872 | 0.3256 |
| Expenditure | 1.2255 | 0.1980 | 6.1901 | 0.0000 |
| ZoneUrban | 63.4600 | 40.0902 | 1.5829 | 0.1162 |
| SexMale | 21.7326 | 15.7808 | 1.3772 | 0.1711 |
| Age2 | 0.0085 | 0.0056 | 1.5297 | 0.1288 |
Por lo anterior, la estimación del valor esperado o predicción queda:
\[ \hat{E}(y_{i}\mid\boldsymbol{x}_{obs,i})=62.2+1.2x_{1i}+63.5x_{2i}+21.7x_{3i}+0.01x_{4i} \] Para calcular la varianza de la estimación, primero se deben obtener las varianzas de la estimación de los parámetros:
vcov(mod_svy)| (Intercept) | Expenditure | ZoneUrban | SexMale | Age2 | |
|---|---|---|---|---|---|
| (Intercept) | 3967.6496 | -11.3101 | 275.9858 | 312.4533 | -0.1852 |
| Expenditure | -11.3101 | 0.0392 | -3.5526 | -0.5004 | 0.0005 |
| ZoneUrban | 275.9858 | -3.5526 | 1607.2280 | -130.7108 | -0.0559 |
| SexMale | 312.4533 | -0.5004 | -130.7108 | 249.0339 | -0.0043 |
| Age2 | -0.1852 | 0.0005 | -0.0559 | -0.0043 | 0.0000 |
Ahora bien, se procede a realizar los cálculos como lo indica la expresión mostrada anteriormente:
xobs <- model.matrix(mod_svy) %>%
data.frame() %>% slice(1) %>% as.matrix()
cov_beta <- vcov(mod_svy) %>% as.matrix()
as.numeric(xobs %*% cov_beta %*% t(xobs))## [1] 1921Si el objetivo ahora es calcular el intervalo de confianza para la predicción se utiliza la siguiente ecuación:
\[ \boldsymbol{x}_{obs,i}\hat{\beta}\pm t_{\left(1-\frac{\alpha}{2},n-p\right)}\sqrt{var\left(\hat{E}\left(y_{i}\mid\boldsymbol{x}_{obs,i}\right)\right)} \]
Para realizar los cálculos en R, se utiliza la función confint y predict como sigue:
pred <- data.frame(predict(mod_svy, type = "link"))
pred_IC <- data.frame(confint(predict(mod_svy, type = "link")))
colnames(pred_IC) <- c("Lim_Inf", "Lim_Sup")
pred_ICAhora, de manera gráfica las predicciones e intervalos se vería de la siguiente manera:
pred <- bind_cols(pred, pred_IC)
pred$Expenditure <- encuesta$Expenditure
pred %>% slice(1:6L)
pd <- position_dodge(width = 0.2)
ggplot(pred %>% slice(1:100L),
aes(x = Expenditure , y = link)) +
geom_errorbar(aes(ymin = Lim_Inf,
ymax = Lim_Sup),
width = .1,
linetype = 1) +
geom_point(size = 2, position = pd) +
theme_bw()Por último, si el interés es hacer una predicción fuera del rango de valores que fue capturado en la muestra. Para esto, supongamos que se desea predecir:
datos_nuevos <- data.frame(Expenditure = 1600,
Age2 = 40^2, Sex = "Male",
Zone = "Urban")La varianza para la predicción se hace siguiendo la siguiente ecuación:
\[ var\left(\hat{E}\left(y_{i}\mid\boldsymbol{x}_{obs,i}\right)\right)=\boldsymbol{x}_{obs,i}^{t}cov\left(\boldsymbol{\beta}\right)\boldsymbol{x}_{obs,i} + \hat{\sigma}^2_{yx} \]
Por tanto, se construye la matriz de observaciones y se calcula la varianza como sigue:
x_noObs = matrix(c(1,1600,1,1,40^2),nrow = 1)
as.numeric(sqrt(x_noObs%*%cov_beta%*%t(x_noObs)))## [1] 244.6Por último, el intervalo de confianza sigue la siguiente ecuación:
\[
\boldsymbol{x}_{obs,i}\hat{\beta}\pm t_{\left(1-\frac{\alpha}{2},n-p\right)}\sqrt{var\left(\hat{E}\left(y_{i}\mid\boldsymbol{x}_{obs,i}\right)\right)+\hat{\sigma}_{yx}^{2}}
\]
En R se hace la predicción de la siguiente manera:
predict(mod_svy, newdata = datos_nuevos, type = "link")## link SE
## 1 2122 245y el intervalo:
confint(predict(mod_svy,newdata = datos_nuevos))| 2.5 % | 97.5 % |
|---|---|
| 1642 | 2601 |