Considerar el conjunto de datos de fecundidad de la ONU (disponible aquí), que contiene la 193 observaciones de tres variables, cada una corresponde a un país. Las variables son:

El objetivo es construir: - Intervalos de confianza para \(\sigma^2\), - Intervalos de confianza individuales y simultános para las componentes de \(\hat{\boldsymbol{\beta}}\), - Regiones de confianza para \(\hat{\boldsymbol{\beta}}\) y sus componentes por pares, - Intervalos de confianza para la respuesta media dado un vector \(\mathbf{x}\), - Intervalos de predicción para una nueva observación.

Cargar los datos en R

file_url <- 'http://users.stat.umn.edu/~sandy/alr3ed/website/data/UN2.txt'
datos <- read.table(file_url, T)
rownames(datos) <- datos$Locality
datos <- datos[,-4]

Graficamos los datos

You must enable Javascript to view this page properly.

Ajuste del modelo

Primero ajustamos el modelo a mano, es decir, con las expresiones que calculamos en clase

\[ \hat{\boldsymbol{\beta}} = \left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\mathbf{y} \qquad \text{y} \qquad \hat{\sigma}^2 = \frac{1}{n-(p+1)}(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})' (\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}}) \]


Declaramos el vector respuesta y la matriz de diseño:

y <- datos[,'logFertility']
X <- cbind(1, datos[,c('logPPgdp', 'Purban')])
X <- as.matrix(X)


Estimamos \(\hat{\boldsymbol{\beta}}\) y \(\sigma^2\):

bh  <- solve(t(X) %*% X) %*% t(X) %*% y; round(bh, 5)
##              [,1]
## 1         2.59300
## logPPgdp -0.12548
## Purban   -0.00352
s2h <- t(y - X%*%bh) %*% (y - X%*%bh) / (193 - 3)
s2h <- as.numeric(s2h); round(s2h, 5)
## [1] 0.15494


Estimamos la varianza de \(\hat{\boldsymbol{\beta}}\):

vbh <- s2h * solve(t(X) %*% X); round(vbh, 6)
##                  1  logPPgdp    Purban
## 1         0.021569 -0.002450  0.000111
## logPPgdp -0.002450  0.000365 -0.000028
## Purban    0.000111 -0.000028  0.000004


Intervalo de confianza para \(\sigma^2\)

Un intervalo de confianza \(100(1-\alpha)\%\) para \(\sigma^2\) está dado por

\[ \left(\frac{(n-p-1)\sigma^2_{MCO}}{\chi^2_{n-p-1}(\alpha/2)}, \frac{(n-p-1)\sigma^2_{MCO}}{\chi^2_{n-p-1}(1-\alpha/2)} \right) \]


Calculamos el intervalo de confianza \(95\%\) para \(\sigma^2\)

gamma1 <- qchisq(0.025, 193 - 3)
gamma2 <- qchisq(0.975, 193 - 3)
s2_LL <- (193 - 3) * s2h / gamma2; s2_LL
## [1] 0.1279598
s2_UL <- (193 - 3) * s2h / gamma1; s2_UL
## [1] 0.1915088

Concluimos que \(\sigma^2\) está contenido en el intervalo \((\) 0.128, 0.192 \()\) con \(95\%\) de confianza.


Intervalos de confianza y pruebas de hipótesis para las componentes de \(\boldsymbol{\beta}\)

Los intervalos de confianza para los componentes de \(\boldsymbol{\beta}\) están dados por

\[ \hat{\beta}_i \pm t_{n - p - 1}(\alpha/2)\sqrt{\hat{V}(\hat{\boldsymbol{\beta}})_{ii}} \]


donde \(\hat{\beta}_i\) es la \(i\)-ésima entrada de \(\hat{\boldsymbol{\beta}}\),
\(t_{n - p - 1}(\alpha/2)\) es el cuantil superior \(\alpha/2\) de una distribución \(t_{n - p - 1}\) y \(\hat{V}(\hat{\boldsymbol{\beta}})_{ii}\) es el \(i\)-ésimo elemento de la diagonal de \(\hat{V}(\hat{\boldsymbol{\beta}})_{ii}\).

Para contrastar las hipótesis \[ H_0: \beta_i = b \qquad \text{vs.} \qquad H_1: \beta_i\neq 0 \]


Se utiliza el estadístico \[ t_i = \frac{\hat{\beta}_i - b}{\sqrt{\hat{V}(\hat{\boldsymbol{\beta}})_{ii}}} \]


La regla de decisión es rechazar si \(\vert{t_i}\vert > t_{n - p - 1}(\alpha/2)\).

Para \(\beta_0\)

b0_LL <- bh[1] - qt(0.975, 93 - 3) * sqrt(vbh[1,1]); b0_LL
## [1] 2.301225
b0_UL <- bh[1] + qt(0.975, 93 - 3) * sqrt(vbh[1,1]); b0_UL
## [1] 2.884766


El intervalo de confianza \(95\%\) para \(\beta_0\) es \((\) 2.301, 2.885 \()\). Podemos estar interesados en decir si el modelo debe o no incluir el intercepto \(\beta_0\), que se puede plantear como un contraste de hipótesis con \(H_0: \beta_0 = 0\). Una forma rápida de hacerlo es notar que el 0 no está incluido en el intervalo de confianza \(95\%\), esto nos sirve para rechazar \(H_0\) con una significancia \(0.05\).

También podemos aplicar el procedimiento usual

t0 <- abs(bh[1]/sqrt(vbh[1,1])); t0 # Estadístico
## [1] 17.65579
tt <- qt(0.975, 93 - 3); tt         # Cuantil de la distribución t  
## [1] 1.986675

Como \(\vert t_0 \vert > t_{n - p - 1}(0.025)\), se rechaza la hipótesis nula con una significancia \(\alpha = 0.05\), es decir, la evidencia sostiene que el modelo incluya a \(\beta_0\).


Para \(\beta_1\)

b1_LL <- bh[2] - qt(0.975, 93 - 3) * sqrt(vbh[2,2]); b1_LL
## [1] -0.1634102
b1_UL <- bh[2] + qt(0.975, 93 - 3) * sqrt(vbh[2,2]); b1_UL
## [1] -0.08754003

El intervalo de confianza \(95\%\) para \(\beta_1\) es \((\) -0.163, -0.088 \()\). Podemos estar interesados en decir si el modelo debe o no incluir el efecto del producto, que se puede plantear como un contraste de hipótesis con \(H_0: \beta_1 = 0\). Nuevamente, una forma rápida de hacerlo es notar que el 0 no está incluido en el intervalo de confianza \(95\%\), esto nos sirve para rechazar \(H_0\) con una significancia \(0.05\). No sólo eso, podemos notar que el intervalo esta conformado por valores negativos, lo que nos lleva sostener que el efecto marginal del producto sobre la fecundidad es negativo, con una significancia \(\alpha = 0.05\).

También podemos aplicar el procedimiento usual

t1 <- abs(bh[2]/sqrt(vbh[2,2])); t1 # Estadístico
## [1] 6.571177
tt <- qt(0.975, 93 - 3); tt         # Cuantil de la distribución t  
## [1] 1.986675

Como \(\vert t_1 \vert > t_{n - p - 1}(0.025)\), se rechaza la hipótesis nula con una significancia \(\alpha = 0.05\), es decir, la evidencia sostiene que el modelo incluya a \(\beta_1\).


Para \(\beta_2\)

b2_LL <- bh[3] - qt(0.975, 93 - 3) * sqrt(vbh[3,3]); b2_LL
## [1] -0.007265906
b2_UL <- bh[3] + qt(0.975, 93 - 3) * sqrt(vbh[3,3]); b2_UL
## [1] 0.0002214697

El intervalo de confianza \(95\%\) para \(\beta_2\) es \((\) -0.00727, 2.210^{-4} \()\). Podemos estar interesados en decir si el modelo debe o no incluir el efecto del porcentaje de población en zonas urbanas, que se puede plantear como un contraste de hipótesis con \(H_0: \beta_2 = 0\). En este caso podemos notar como el 0 está incluido en el intervalo de confianza \(95\%\), esto nos indica que \(H_0\) no se rechaza con una significancia \(0.05\). Esto es evidencia para sostener que el porcentaje de población en zonas urbans no tiene efecto sobre la fecundidad, con una significancia \(\alpha = 0.05\).

También podemos aplicar el procedimiento usual

t2 <- abs(bh[3]/sqrt(vbh[3,3])); t2 # Estadístico
## [1] 1.869147
tt <- qt(0.975, 93 - 3); tt         # Cuantil de la distribución t  
## [1] 1.986675

Como \(\vert t_2 \vert < t_{n - p - 1}(0.025)\), no se rechaza la hipótesis nula con una significancia \(\alpha = 0.05\), es decir, la evidencia sostiene que el modelo no debería incluir a \(\beta_2\).

Con esto terminamos las pruebas de hipótesis e intervalos de confianza individuales sobre las componentes de \(\boldsymbol{\beta}\), después de hacer algunos comentarios pocederemos a hacer pruebas e intervalos simultáneos.


Gráficas de regresión parcial

El modelo simple de logFertility sobre Purban indica que esta variable sí tiene un efecto sobre la respuesta. Sin embargo, cuando incluimos la variable logPPgdp el efecto de Purban ya no es significativo. Esto se debe a que las variables logPPgdd y Purban están relacionadas, por lo que casi toda la variabilidad de logFertility que es explicada por Purban está contenida en la variabilidad explicada por logPPgdp. ¿Cómo podemos ver esto?

Primero consideremos los gráficos de dispersión por pares de variables.


Vemos parece que ambas variables tienen una buena relación con logFertility, pero están relacionadas entre ellas. Procedemos de forma secuencial a ajustar el modelo. Primero ajustamos logFertility contra logPPgdp.

modelo1 <- lm(logFertility ~ logPPgdp, datos)
summary(modelo1)
## 
## Call:
## lm(formula = logFertility ~ logPPgdp, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11877 -0.18762  0.07041  0.26082  0.90101 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.70322    0.13538   19.97   <2e-16 ***
## logPPgdp    -0.15330    0.01204  -12.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3962 on 191 degrees of freedom
## Multiple R-squared:  0.4591, Adjusted R-squared:  0.4563 
## F-statistic: 162.1 on 1 and 191 DF,  p-value: < 2.2e-16


Los residuos del modelo anterior contienen la variabilidad de logFertility no explicada por logPPgdp. Al igual que si ajustamos un modelo de Purban contra logPPgdp, los residuales contentran la variabilidad en Purban que no es explicada por logPPgdp.

modelo2 <- lm(Purban ~ logPPgdp, datos)
summary(modelo2)
## 
## Call:
## lm(formula = Purban ~ logPPgdp, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.777  -8.692   0.635  10.150  38.851 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -31.2934     5.1647  -6.059 7.15e-09 ***
## logPPgdp      7.8988     0.4593  17.198  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.11 on 191 degrees of freedom
## Multiple R-squared:  0.6076, Adjusted R-squared:  0.6056 
## F-statistic: 295.8 on 1 and 191 DF,  p-value: < 2.2e-16


Finalmente, ajustamos un modelo para explicar la variabilidad residual de logFertility a partir de la varianza residual de Purban, es decir, cuando en ambas variables eliminamos el efecto de logPPgdp.

res_logfert <- residuals(modelo1)
res_Purban <- residuals(modelo2)
modelo3 <- lm(res_logfert ~ res_Purban)
summary(modelo3)
## 
## Call:
## lm(formula = res_logfert ~ res_Purban)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05378 -0.16952  0.06835  0.25587  0.89290 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -1.575e-18  2.826e-02   0.000   1.0000  
## res_Purban  -3.522e-03  1.879e-03  -1.874   0.0624 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3926 on 191 degrees of freedom
## Multiple R-squared:  0.01806,    Adjusted R-squared:  0.01291 
## F-statistic: 3.512 on 1 and 191 DF,  p-value: 0.06245


El intercepto de este modelo es 0, ya que ambas variables son residuales de modelos con intercepto, por lo que su promedio es 0. El coeficiente de este modelo, es el coeficiente que se obtiene para Purban en el modelo que incluye las dos variables explicativas. Lo anterior funciona si cambiamos logPPgdp por Purban y viceversa.

La gráfica de regresión parcial es la que se obtiene de gráficar los residuales ajustados en el modelo 3.

plot(res_Purban, res_logfert, pch = 16, col = 'steelblue',
     xlab = 'Purban | logPPgdp', ylab = 'logFertility | logPPgdp')


A partir de la gráfica anterior podemos entender de una mejor manera cuál es la relación de la variable respuesta y cada variable explicativa, en presencia de las demás variables explicativas. Podemos identificar desviaciones al supuesto de normalidad, a la especificación del modelo, observaciones influyentes, entre otras. Las gráficas de variable añadida se pueden construir para cada variable de una regresión múltiple con cualquier número de variables.


Análisis de varianza

Ahora procedemos a hacer el análisis de varianza de nuestro modelo de regresión. La parte de la tabla ANOVA que podemos llenar sin mayor complicación es


F.V. G.L. S.C. C.M. F
Regresión 2
Error 190
Total corregido 192


La suma de cuadrados del total corregido se obtiene directamente de las observaciones. La suma de cuadrados del error la obtenemos del a partir de los valores ajustados. La suma de cuadrados de regresión la obtenemos por diferencia.

sc_tc <- sum((y - mean(y))^2);    sc_tc
## [1] 55.43156
sc_error <- sum((y - X%*%bh)^2);  sc_error
## [1] 29.439
sc_reg <- sc_tc - sc_error;       sc_reg
## [1] 25.99257


El cuadrado medio (C.M.) se obtiene dividiendo las sumas de cuadrados entre los grados de libertad que corresponde. Finalmente el estadístico \(F\) se obtiene de dividir \(CM_{reg}\) entre \(CM_{error}\).


F.V. G.L. S.C. C.M. F
Regresión 2 25.992 12.996 83.899
Error 190 29.439 0.154
Total corregido 192 55.431


El cuantil superior \(0.05\) de una distribución \(F\) con 2 y 190 grados de libertad es 3.04, por lo que hay evidencia en contra de la hipótesis nula, es decir, al menos una de las variables del modelo tiene algún efecto en la fecundidad.


Región de confianza para \(\boldsymbol{\beta}\)

Una región de confianza \(100(1-\alpha)\%\) para \(\boldsymbol{\beta}\) está dada por

\[ \left\lbrace \boldsymbol{\beta} \in \mathbb{R}^{p+1} \mid (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) '\mathbf{X}'\mathbf{X} (\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \leq \hat{\sigma}^2_{MCO}(p+1)K_\alpha \right\rbrace \]


con \(K_\alpha = F_{p+1, n-p-1}(\alpha)\).

En este caso utilizaremos la función ellipse3d para únicamente graficar la región de confianza. Más adelante veremos un ejemplo con regresión simple en donde calcularemos con detalle la región de confianza para los parámetros.


You must enable Javascript to view this page properly.


Región de confianza para \(\beta_1\) y \(\beta_2\)

Para calcular regiones de confianza para \(\beta_1\) y \(\beta_2\) utilizamos la distribución marginal de \(\hat{\beta}_1\) y \(\hat{\beta}_2\).

La siguiente función corresponde a la forma cuadrátrica que define la región de confianza.

region2 <- function(x, y, bh, X){
  b0 <- matrix(c(x, y), 2, 1)
  b <- bh[2:3]
  aux <- solve(t(X)%*%X)
  m <- solve(aux[2:3, 2:3])
  res <- t(b - b0)%*%m%*%(b-b0)
  return(res)
}


Se crea una malla de puntos en lo que se evalua la función anterior.

seqx <- seq(bh[2] - 0.06, bh[2] + 0.06, length.out = 250)
seqy <- seq(bh[3]-0.006, bh[3]+0.006, length.out = 250)

m_aux <- matrix(0, 250, 250)
for (i in 1:250){
  for (j in 1:250){
    m_aux[i, j] <- region2(seqx[i], seqy[j], bh, X)
  }
}


Por último se grafican las regiones de confianza 90%, 95% y 99%.

confianza <- c(2 * s2h * qf(0.90, 2, 190),
               2 * s2h * qf(0.95, 2, 190),
               2 * s2h * qf(0.99, 2, 190))
image(seqx, seqy, m_aux, breaks = c(0, confianza), xlim = c(-0.2, 0),
      col = c('steelblue3', 'steelblue2', 'steelblue1'),
      xlab = 'beta1', ylab = 'beta2')
contour(seqx, seqy, m_aux, levels = confianza, drawlabels = F, add = T)
legend('topright', legend = c('90%', '95%', '99%'), pch = 16,
       col = c('steelblue3', 'steelblue2', 'steelblue1'), cex = 1.5)

Las regiones de confianza anteriores sirven para contrastar la hipótesis que se plantea en el ANOVA. Si la región de confianza no contiene al punto \((0, 0)\), la prueba ANOVA se rechaza con el nivel de significancia específicado.


Intervalos de confianza simultáneos

Ahora construiremos intervalos de confianza simultáneos para las componentes de \(\boldsymbol{\beta}\). Estos intervalos nos permiten hacer pruebas de hipótesis simúltáneas para los coeficientes del modelo, a diferencia de las pruebas \(t\) que son individuales. El propósito también es comparar la longitud de los intervalos.


Bonferroni

El método de Bonferroni para intervalos secundarios se basa en el principio de inclusión y exclusión, que se traduce en cambiar los cuantiles superiores \(\alpha/2\) por \(\alpha/2k\), donde \(k\) es el número de intervalos que se quieren construir, en nuestro caso \(k=3\).

\[ \hat{\beta}_i \pm t_{n - p - 1}(\alpha/2k)\sqrt{\hat{V}(\hat{\boldsymbol{\beta}})_{ii}} \]


Para \(\beta_0\)

b0_LL_B <- bh[1] - qt(1 - 0.05/6, 193 - 3) * sqrt(vbh[1,1]); b0_LL_B
## [1] 2.238267
b0_UL_B <- bh[1] + qt(1 - 0.05/6, 193 - 3) * sqrt(vbh[1,1]); b0_UL_B
## [1] 2.947725


Para \(\beta_1\)

b1_LL_B <- bh[2] - qt(1 - 0.05/6, 193 - 3) * sqrt(vbh[2,2]); b1_LL_B
## [1] -0.1715959
b1_UL_B <- bh[2] + qt(1 - 0.05/6, 193 - 3) * sqrt(vbh[2,2]); b1_UL_B
## [1] -0.07935432


Para \(\beta_2\)

b2_LL_B <- bh[3] - qt(1 - 0.05/6, 193 - 3) * sqrt(vbh[3,3]); b2_LL_B
## [1] -0.008073726
b2_UL_B <- bh[3] + qt(1 - 0.05/6, 193 - 3) * sqrt(vbh[3,3]); b2_UL_B
## [1] 0.001029289


Hotelling-Scheffé

El método de Hotelling-Scheffé se basa en maximizar la forma cuadrática del estadístico utilizado en los intervalos individuales. Esto se traduce en cambiar el cuantil de la \(t\) por un múltiplo del cuantil de la \(F\). La ventaja de este método es que es válido para cualquier número de intervalos simultáneos.

\[ \hat{\beta}_i \pm \sqrt{K_\alpha \hat{V}(\hat{\boldsymbol{\beta}})_{ii}} \]

con \(K_\alpha = (p+1)F_{p+1, n-p-1}(\alpha)\).


Para \(\beta_0\)

b0_LL_HS <- bh[1] - sqrt(3 * qf(0.95, 3, 190)) * sqrt(vbh[1,1]); b0_LL_HS
## [1] 2.178736
b0_UL_HS <- bh[1] + sqrt(3 * qf(0.95, 3, 190)) * sqrt(vbh[1,1]); b0_UL_HS
## [1] 3.007256


Para \(\beta_1\)

b1_LL_HS <- bh[2] - sqrt(3 * qf(0.95, 3, 190)) * sqrt(vbh[2,2]); b1_LL_HS
## [1] -0.179336
b1_UL_HS <- bh[2] + sqrt(3 * qf(0.95, 3, 190)) * sqrt(vbh[2,2]); b1_UL_HS
## [1] -0.07161428


Para \(\beta_2\)

b2_LL_HS <- bh[3] - sqrt(3 * qf(0.95, 3, 190)) * sqrt(vbh[3,3]); b2_LL_HS
## [1] -0.008837565
b2_UL_HS <- bh[3] + sqrt(3 * qf(0.95, 3, 190)) * sqrt(vbh[3,3]); b2_UL_HS
## [1] 0.001793128


Comaparación de los intervalos de confianza


Región e intervalos simultáneos de confianza para \(\beta_1\) y \(\beta_2\)

¿Qué relación tiene la región de confianza con los intervalos simultáneos? Ambos definen conjuntos de valores plausibles para el vector de parámetros desconocido. La región de confianza se basa en la distribución conjunta del vector de estimadores. Tiene la ventaja de tomar en cuenta la estructura de convarianzas, por lo que el conjunto de valores plausibles que define tiene el menor tamaño posible. Tiene la desventaja que para dimensiones mayores a 3 es difícil de visualizar.

Los intervalos simultáneos se basan en algún método ajuste de los intervalos individuales. El método de Bonferroni utiliza el principio de inclusión-exclusión. El método de Hotelling-Scheffé se basa en un argumento de optimización. Tienen la ventaja que definen regiones fáciles de interpretar (\(n\)-cajas). Tienen la desventaja que definen regiones plausibles muy amplias. Como los métodos de ajuste no son exactos resultan muy conservadores en algunos casos, lo que puede llevar a contradicciones si se comparar contra la región de confianza.

Para ilustrar lo anterior consideramos únicamente los parámetros \(\beta_1\) y \(\beta_2\). Ya habíamos calculado una región de confianza para el vector \((\beta_1, \beta_2)\), sólo hace falta calcular intervalos simultáneos únicamente para estos dos parámetros.


Bonferroni

Utililizamos la formula anterior, pero ahora con \(k=2\):

b1_LL_B2 <- bh[2] - qt(1 - 0.05/4, 193 - 3) * sqrt(vbh[2,2]); b1_LL_B2
## [1] -0.168616
b1_UL_B2 <- bh[2] + qt(1 - 0.05/4, 193 - 3) * sqrt(vbh[2,2]); b1_UL_B2
## [1] -0.08233421
b2_LL_B2 <- bh[3] - qt(1 - 0.05/4, 193 - 3) * sqrt(vbh[3,3]); b2_LL_B2
## [1] -0.00777965
b2_UL_B2 <- bh[3] + qt(1 - 0.05/4, 193 - 3) * sqrt(vbh[3,3]); b2_UL_B2
## [1] 0.000735214


Hotelling-Scheffé

Ajustamos el valor de \(K_\alpha\) ya que \(\beta_0\) no es de nuestro interés para esta comparación:

b1_LL_HS2 <- bh[2] - sqrt(2 * qf(0.95, 2, 190)) * sqrt(vbh[2,2])
b1_UL_HS2 <- bh[2] + sqrt(2 * qf(0.95, 2, 190)) * sqrt(vbh[2,2])
b2_LL_HS2 <- bh[3] - sqrt(2 * qf(0.95, 2, 190)) * sqrt(vbh[3,3])
b2_UL_HS2 <- bh[3] + sqrt(2 * qf(0.95, 2, 190)) * sqrt(vbh[3,3])


En la siguiente gráfica se muestran los intervalos simultáneos y la región de confianza.


Ajuste del modelo

Calculamos el coeficiente \(R^2\)

r2 <- 1 - sc_error/sc_tc; r2
## [1] 0.4689128


Calculamos el \(R^2\) ajustado

r2adj <- 1 - (sc_error / 190) / (sc_tc / 192); r2adj
## [1] 0.4633224


Comparamos los resultados con la salida de R

ajuste1 <- lm(logFertility ~ logPPgdp + Purban, datos)
summary(ajuste1)
## 
## Call:
## lm(formula = logFertility ~ logPPgdp + Purban, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05378 -0.16952  0.06835  0.25587  0.89290 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.592996   0.146864  17.656  < 2e-16 ***
## logPPgdp    -0.125475   0.019095  -6.571 4.67e-10 ***
## Purban      -0.003522   0.001884  -1.869   0.0631 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3936 on 190 degrees of freedom
## Multiple R-squared:  0.4689, Adjusted R-squared:  0.4633 
## F-statistic: 83.88 on 2 and 190 DF,  p-value: < 2.2e-16


Los resultados coinciden.