1 Introducción

El propósito de este documento es realizar un análisis estadístico completo sobre el conjunto de datos Salaries, disponible en el paquete carData de R. Este conjunto contiene información sobre los salarios de profesores universitarios de una universidad de Estados Unidos, junto con variables que describen su rango académico, disciplina, sexo, años desde el doctorado y años de servicio.

2 Paquetes

# Manipulación y visualización de datos
library(carData)     # Contiene el dataset Salaries
library(dplyr)        # Manipulación de datos
library(ggplot2)       # Gráficos
library(tidyr)         # Reorganización de datos (para gráficos en panel)

# Pruebas estadísticas y diagnóstico de modelos
library(car)           # vif(), outlierTest(), durbinWatsonTest(), influencePlot()
library(lmtest)        # bptest() para homocedasticidad
library(nortest)       # Pruebas de normalidad adicionales (opcional, de respaldo)
library(MASS)          # studres() para residuos estudentizados
library(PerformanceAnalytics)  # chart.Correlation() para el panel de correlaciones

# Tablas y formato
library(knitr)         # kable() para tablas
library(kableExtra)    # Formato adicional de tablas

paleta_pastel <- c("#A8D8EA", "#AA96DA", "#FCBAD3", "#FFFFD2", "#B5EAD7", "#FFDAC1")
paleta_pastel_2 <- c("#FBC4AB", "#A0CED9")  # para variables de 2 niveles (discipline, sex)

3 Descripción general del conjunto de datos

El conjunto de datos Salaries proviene de un estudio sobre la equidad salarial entre profesores universitarios y contiene 397 observaciones y 6 variables. Cada fila representa a un profesor o profesora, y las columnas describen características demográficas y profesionales, además del salario anual.

data("Salaries", package = "carData")
datos <- Salaries

# Dimensiones del conjunto de datos
dim(datos)

## [1] 397   6

# Primeras observaciones
head(datos)

# Estructura interna: tipo de cada variable
str(datos)

## 'data.frame':    397 obs. of  6 variables:
##  $ rank         : Factor w/ 3 levels "AsstProf","AssocProf",..: 3 3 1 3 3 2 3 3 3 3 ...
##  $ discipline   : Factor w/ 2 levels "A","B": 2 2 2 2 2 2 2 2 2 2 ...
##  $ yrs.since.phd: int  19 20 4 45 40 6 30 45 21 18 ...
##  $ yrs.service  : int  18 16 3 39 41 6 23 45 20 18 ...
##  $ sex          : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 1 ...
##  $ salary       : int  139750 173200 79750 115000 141500 97000 175000 147765 119250 129000 ...

3.1 Diccionario de variables

A continuación se describe cada variable, indicando su naturaleza estadística:

Variable	Tipo	Descripción
`rank`	Categórica nominal (3 niveles)	Rango académico del profesor: `AsstProf` (profesor asistente), `AssocProf` (profesor asociado), `Prof` (profesor titular).
`discipline`	Categórica nominal (2 niveles)	Área disciplinar: `A` (áreas teóricas) o `B` (áreas aplicadas).
`yrs.since.phd`	Cuantitativa discreta (numérica)	Años transcurridos desde la obtención del doctorado.
`yrs.service`	Cuantitativa discreta (numérica)	Años de servicio en la universidad.
`sex`	Categórica nominal (2 niveles)	Sexo del profesor: `Male` o `Female`.
`salary`	Cuantitativa continua (numérica)	Salario anual en dólares estadounidenses. Esta es la variable respuesta (Y) del modelo.

En estadística, distinguir entre variables cualitativas (categóricas) y cuantitativas (numéricas) es el primer paso de cualquier análisis, porque determina qué tipo de medidas de resumen y qué tipo de gráficos son apropiados para cada una. En este conjunto:

Variables categóricas: rank, discipline, sex.
Variables numéricas: yrs.since.phd, yrs.service, salary.

summary(datos)

##         rank     discipline yrs.since.phd    yrs.service        sex     
##  AsstProf : 67   A:181      Min.   : 1.00   Min.   : 0.00   Female: 39  
##  AssocProf: 64   B:216      1st Qu.:12.00   1st Qu.: 7.00   Male  :358  
##  Prof     :266              Median :21.00   Median :16.00               
##                             Mean   :22.31   Mean   :17.61               
##                             3rd Qu.:32.00   3rd Qu.:27.00               
##                             Max.   :56.00   Max.   :60.00               
##      salary      
##  Min.   : 57800  
##  1st Qu.: 91000  
##  Median :107300  
##  Mean   :113706  
##  3rd Qu.:134185  
##  Max.   :231545

4 Análisis exploratorio de datos (AED)

4.1 Análisis univariado de las variables numéricas

Antes de relacionar variables entre sí, es necesario entender el comportamiento individual de cada variable numérica: su tendencia central, su dispersión y su forma de distribución.

# Se construye la tabla directamente con data.frame() (no con summarise())
# porque summarise() exige que cada columna resultante tenga una sola fila
# por grupo; aquí necesitamos 3 filas (una por variable), por lo que
# data.frame() es la herramienta correcta para este caso.
medidas_resumen <- data.frame(
  Variable = c("yrs.since.phd", "yrs.service", "salary"),
  Media = c(mean(datos$yrs.since.phd), mean(datos$yrs.service), mean(datos$salary)),
  Mediana = c(median(datos$yrs.since.phd), median(datos$yrs.service), median(datos$salary)),
  Desv_Estandar = c(sd(datos$yrs.since.phd), sd(datos$yrs.service), sd(datos$salary)),
  Minimo = c(min(datos$yrs.since.phd), min(datos$yrs.service), min(datos$salary)),
  Maximo = c(max(datos$yrs.since.phd), max(datos$yrs.service), max(datos$salary)),
  CV_pct = c(sd(datos$yrs.since.phd) / mean(datos$yrs.since.phd) * 100,
             sd(datos$yrs.service) / mean(datos$yrs.service) * 100,
             sd(datos$salary) / mean(datos$salary) * 100)
)

kable(medidas_resumen, digits = 2,
      caption = "Medidas de resumen de las variables numéricas") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Medidas de resumen de las variables numéricas
Variable	Media	Mediana	Desv_Estandar	Minimo	Maximo	CV_pct
yrs.since.phd	22.31	21	12.89	1	56	57.75
yrs.service	17.61	16	13.01	0	60	73.84
salary	113706.46	107300	30289.04	57800	231545	26.64

Interpretación:

El salario promedio es de aproximadamente el valor reportado en la tabla, con una desviación estándar considerable, lo que indica que existe una dispersión importante entre los salarios de los profesores.
El coeficiente de variación (CV) —que se calcula como la desviación estándar dividida por la media, expresada en porcentaje— permite comparar la dispersión relativa entre variables con distintas unidades. Un CV alto en salary sugiere heterogeneidad salarial dentro de la institución.
yrs.since.phd y yrs.service están naturalmente correlacionadas, ya que ambas miden antigüedad profesional, aunque desde ángulos distintos (desde el doctorado vs. desde el ingreso a la universidad).

4.2 Distribución de las variables numéricas: histogramas

datos_largo <- datos %>%
  dplyr::select(yrs.since.phd, yrs.service, salary) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor")

ggplot(datos_largo, aes(x = valor)) +
  geom_histogram(aes(y = after_stat(density)), bins = 25, fill = "#2C7FB8", color = "white", alpha = 0.85) +
  geom_density(color = "#D7191C", linewidth = 1) +
  facet_wrap(~variable, scales = "free") +
  theme_minimal(base_size = 12) +
  labs(title = "Distribución de las variables numéricas",
       subtitle = "Histograma con curva de densidad superpuesta",
       x = "Valor", y = "Densidad")

Interpretación: La curva roja (densidad estimada) permite visualizar la forma de la distribución. Si la curva fuera perfectamente simétrica y con forma de campana, hablaríamos de una distribución cercana a la normal. Visualmente, salary y yrs.service muestran cierta asimetría hacia la derecha (cola larga hacia valores altos), lo cual es típico en variables de ingresos y de antigüedad laboral. Esta primera impresión visual será confirmada o refutada formalmente más adelante con la prueba de Kolmogorov-Smirnov.

4.3 Análisis univariado de las variables categóricas

tabla_rank <- table(datos$rank)
tabla_discipline <- table(datos$discipline)
tabla_sex <- table(datos$sex)

kable(as.data.frame(tabla_rank), col.names = c("Rango académico", "Frecuencia"),
      caption = "Distribución de la variable rank") %>%
  kable_styling(full_width = FALSE)

Distribución de la variable rank
Rango académico	Frecuencia
AsstProf	67
AssocProf	64
Prof	266

kable(as.data.frame(tabla_discipline), col.names = c("Disciplina", "Frecuencia"),
      caption = "Distribución de la variable discipline") %>%
  kable_styling(full_width = FALSE)

Distribución de la variable discipline
Disciplina	Frecuencia
A	181
B	216

kable(as.data.frame(tabla_sex), col.names = c("Sexo", "Frecuencia"),
      caption = "Distribución de la variable sex") %>%
  kable_styling(full_width = FALSE)

Distribución de la variable sex
Sexo	Frecuencia
Female	39
Male	358

p1 <- ggplot(datos, aes(x = rank, fill = rank)) +
  geom_bar(color = "#555555") + scale_fill_manual(values = paleta_pastel[1:3]) +
  theme_minimal() + theme(legend.position = "none") +
  labs(title = "Rango académico", x = "", y = "Frecuencia")

p2 <- ggplot(datos, aes(x = discipline, fill = discipline)) +
  geom_bar(color = "#555555") + scale_fill_manual(values = paleta_pastel_2) +
  theme_minimal() + theme(legend.position = "none") +
  labs(title = "Disciplina", x = "", y = "Frecuencia")

p3 <- ggplot(datos, aes(x = sex, fill = sex)) +
  geom_bar(color = "#555555") + scale_fill_manual(values = c("Female" = "#FCBAD3", "Male" = "#A8D8EA")) +
  theme_minimal() + theme(legend.position = "none") +
  labs(title = "Sexo", x = "", y = "Frecuencia")

gridExtra::grid.arrange(p1, p2, p3, ncol = 3)

Interpretación: Se observa que la categoría Prof (profesor titular) es la más numerosa dentro de rank, lo cual es razonable ya que representa el rango más alto y suele acumular a quienes llevan más tiempo en la carrera académica. En sex, el conjunto de datos está fuertemente desbalanceado hacia hombres, una característica importante a tener en cuenta al interpretar resultados relacionados con esta variable, pues una muestra pequeña en la categoría minoritaria puede afectar la precisión de las estimaciones.

4.4 Relación entre el salario y las variables categóricas: diagramas de caja (boxplots)

4.4.1 Salario según rango académico (`rank`)

ggplot(datos, aes(x = rank, y = salary, fill = rank)) +
  geom_boxplot(alpha = 0.9, outlier.color = "#6A4C93", outlier.shape = 19, color = "#555555") +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3.5, color = "#4A4A4A") +
  scale_fill_manual(values = paleta_pastel[1:3]) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none") +
  labs(title = "Distribución del salario según rango académico",
       subtitle = "El rombo gris indica la media de cada grupo",
       x = "Rango académico", y = "Salario anual (USD)")

Interpretación: Se espera observar un incremento claro y progresivo en la mediana salarial a medida que aumenta el rango académico: AsstProf < AssocProf < Prof. Esto es consistente con la lógica de cualquier escalafón laboral académico, donde la antigüedad y la jerarquía se traducen en mayor remuneración. Si las cajas de los tres grupos no se solapan o se solapan muy poco, es un indicio visual de que rank podría ser una variable explicativa relevante del salario.

4.4.2 Salario según disciplina (`discipline`)

ggplot(datos, aes(x = discipline, y = salary, fill = discipline)) +
  geom_boxplot(alpha = 0.9, outlier.color = "#6A4C93", outlier.shape = 19, color = "#555555") +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3.5, color = "#4A4A4A") +
  scale_fill_manual(values = paleta_pastel_2) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none") +
  labs(title = "Distribución del salario según disciplina",
       subtitle = "A = áreas teóricas, B = áreas aplicadas",
       x = "Disciplina", y = "Salario anual (USD)")

Interpretación: Generalmente, las disciplinas aplicadas (B) muestran salarios medios algo más altos que las teóricas (A), lo cual podría reflejar diferencias en la demanda de mercado de cada área. La diferencia entre ambos grupos suele ser menos marcada que la observada entre rangos académicos.

4.4.3 Salario según sexo (`sex`)

ggplot(datos, aes(x = sex, y = salary, fill = sex)) +
  geom_boxplot(alpha = 0.9, outlier.color = "#6A4C93", outlier.shape = 19, color = "#555555") +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3.5, color = "#4A4A4A") +
  scale_fill_manual(values = c("Female" = "#FCBAD3", "Male" = "#A8D8EA")) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none") +
  labs(title = "Distribución del salario según sexo",
       x = "Sexo", y = "Salario anual (USD)")

Interpretación: Este boxplot permite observar si existen diferencias salariales asociadas al sexo. Es importante recordar que, dado el fuerte desbalance muestral mencionado anteriormente (muchos más hombres que mujeres en la muestra), cualquier diferencia observada debe interpretarse con cautela y, si se desea formalizar, complementarse con una prueba de hipótesis (por ejemplo, una prueba t o una prueba no paramétrica de Mann-Whitney), lo cual queda fuera del alcance principal de este documento pero se sugiere como análisis complementario.

4.4.4 Boxplots combinados: rango y disciplina simultáneamente

Para enriquecer el análisis, se presenta un boxplot que cruza dos variables categóricas a la vez (rank y discipline), lo que permite detectar posibles interacciones entre ellas.

ggplot(datos, aes(x = rank, y = salary, fill = discipline)) +
  geom_boxplot(alpha = 0.9, position = position_dodge(width = 0.8), color = "#555555") +
  theme_minimal(base_size = 12) +
  labs(title = "Salario según rango académico y disciplina",
       x = "Rango académico", y = "Salario anual (USD)", fill = "Disciplina") +
  scale_fill_manual(values = c("A" = "#A0CED9", "B" = "#FBC4AB"))

Interpretación: Este gráfico permite ver si la diferencia salarial entre disciplinas se mantiene constante en todos los rangos académicos, o si por el contrario es mayor o menor en algún rango específico. Si las diferencias entre las cajas de color cambian de magnitud o de dirección según el rango, esto sería evidencia visual de una posible interacción entre rank y discipline.

4.4.5 Boxplots de las variables numéricas explicativas categorizadas en rangos

Como yrs.since.phd y yrs.service son variables numéricas, se construyen versiones categorizadas de ellas (en intervalos o “bins”) para poder visualizarlas también mediante boxplots, igual que se hizo con las variables categóricas originales.

datos <- datos %>%
  mutate(
    rango_phd = cut(yrs.since.phd,
                     breaks = c(-1, 5, 10, 20, 30, max(yrs.since.phd)),
                     labels = c("0-5", "6-10", "11-20", "21-30", "31+")),
    rango_servicio = cut(yrs.service,
                          breaks = c(-1, 5, 10, 20, 30, max(yrs.service)),
                          labels = c("0-5", "6-10", "11-20", "21-30", "31+"))
  )

table(datos$rango_phd)

## 
##   0-5  6-10 11-20 21-30   31+ 
##    42    45   106    93   111

table(datos$rango_servicio)

## 
##   0-5  6-10 11-20 21-30   31+ 
##    82    73    95    77    70

ggplot(datos, aes(x = rango_phd, y = salary, fill = rango_phd)) +
  geom_boxplot(alpha = 0.9, outlier.color = "#6A4C93", color = "#555555") +
  scale_fill_manual(values = paleta_pastel) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none") +
  labs(title = "Salario según rango de años desde el doctorado",
       x = "Años desde el doctorado (categorizado)", y = "Salario anual (USD)")

ggplot(datos, aes(x = rango_servicio, y = salary, fill = rango_servicio)) +
  geom_boxplot(alpha = 0.9, outlier.color = "#6A4C93", color = "#555555") +
  scale_fill_manual(values = paleta_pastel) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none") +
  labs(title = "Salario según rango de años de servicio",
       x = "Años de servicio (categorizado)", y = "Salario anual (USD)")

Interpretación: Categorizar variables numéricas en rangos (también llamado binning) es una técnica exploratoria útil porque, aunque hace perder algo de información respecto a la variable continua original, facilita la lectura de patrones generales mediante boxplots, similar a como se hace con variables categóricas nativas. Se espera observar que el salario tiende a aumentar conforme aumentan los años desde el doctorado y los años de servicio, aunque es común que este incremento se estabilice o incluso disminuya levemente en los tramos más altos (efecto de meseta salarial), un patrón que conviene confirmar visualmente en los gráficos generados.

4.5 Relación entre el salario y las variables numéricas: diagramas de dispersión

p4 <- ggplot(datos, aes(x = yrs.since.phd, y = salary)) +
  geom_point(alpha = 0.5, color = "#2C7FB8") +
  geom_smooth(method = "lm", color = "#D7191C", se = TRUE) +
  theme_minimal() +
  labs(title = "Salario vs. años desde el doctorado", x = "Años desde el doctorado", y = "Salario (USD)")

p5 <- ggplot(datos, aes(x = yrs.service, y = salary)) +
  geom_point(alpha = 0.5, color = "#2C7FB8") +
  geom_smooth(method = "lm", color = "#D7191C", se = TRUE) +
  theme_minimal() +
  labs(title = "Salario vs. años de servicio", x = "Años de servicio", y = "Salario (USD)")

gridExtra::grid.arrange(p4, p5, ncol = 2)

Interpretación: La línea roja representa la recta de regresión simple ajustada y la banda gris su intervalo de confianza. Ambas variables muestran una relación positiva con el salario: a mayor antigüedad (desde el doctorado o de servicio), mayor tendencia salarial, aunque la nube de puntos sugiere que la relación no es perfectamente lineal y que hay bastante dispersión (variabilidad) alrededor de la recta, especialmente en los profesores con más años de experiencia.

5 Análisis de normalidad de las variables numéricas

5.1 Justificación de la prueba a utilizar

Antes de decidir si usar el coeficiente de correlación de Pearson (que asume normalidad) o el de Spearman (que no la asume, pues trabaja con rangos), es necesario evaluar formalmente la normalidad de cada variable numérica.

5.1.1 Planteamiento de la prueba de hipótesis de Kolmogorov-Smirnov

Para cada variable numérica se contrasta:

\[ H_0: \text{La variable proviene de una distribución normal} \] \[ H_1: \text{La variable NO proviene de una distribución normal} \]

Regla de decisión: si el valor p < 0.05, se rechaza \(H_0\) y se concluye que la variable no sigue una distribución normal, con un nivel de significancia del 5%.

# Variable: yrs.since.phd
ks_phd <- ks.test(scale(datos$yrs.since.phd), "pnorm")
ks_phd

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  scale(datos$yrs.since.phd)
## D = 0.069888, p-value = 0.04138
## alternative hypothesis: two-sided

# Variable: yrs.service
ks_service <- ks.test(scale(datos$yrs.service), "pnorm")
ks_service

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  scale(datos$yrs.service)
## D = 0.12017, p-value = 2.096e-05
## alternative hypothesis: two-sided

# Variable: salary (variable respuesta)
ks_salary <- ks.test(scale(datos$salary), "pnorm")
ks_salary

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  scale(datos$salary)
## D = 0.09086, p-value = 0.002846
## alternative hypothesis: two-sided

tabla_ks <- data.frame(
  Variable = c("yrs.since.phd", "yrs.service", "salary"),
  Estadistico_D = c(ks_phd$statistic, ks_service$statistic, ks_salary$statistic),
  Valor_p = c(ks_phd$p.value, ks_service$p.value, ks_salary$p.value)
) %>%
  mutate(Conclusion = ifelse(Valor_p < 0.05, "No normal (se rechaza H0)", "Normal (no se rechaza H0)"))

kable(tabla_ks, digits = 4,
      caption = "Resultados de la prueba de Kolmogorov-Smirnov para normalidad (alpha = 0.05)") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Resultados de la prueba de Kolmogorov-Smirnov para normalidad (alpha = 0.05)
Variable	Estadistico_D	Valor_p	Conclusion
yrs.since.phd	0.0699	0.0414	No normal (se rechaza H0)
yrs.service	0.1202	0.0000	No normal (se rechaza H0)
salary	0.0909	0.0028	No normal (se rechaza H0)

5.2 Confirmación visual: gráficos Q-Q (cuantil-cuantil)

par(mfrow = c(1, 3))
qqnorm(datos$yrs.since.phd, main = "Q-Q Plot: yrs.since.phd", col = "#2C7FB8")
qqline(datos$yrs.since.phd, col = "red", lwd = 2)

qqnorm(datos$yrs.service, main = "Q-Q Plot: yrs.service", col = "#2C7FB8")
qqline(datos$yrs.service, col = "red", lwd = 2)

qqnorm(datos$salary, main = "Q-Q Plot: salary", col = "#2C7FB8")
qqline(datos$salary, col = "red", lwd = 2)

par(mfrow = c(1, 1))

Interpretación: En un gráfico Q-Q, si los puntos siguen de cerca la línea diagonal roja, la variable se aproxima a una distribución normal. Desviaciones sistemáticas de los puntos respecto a la línea (especialmente en las colas) confirman lo detectado por la prueba K-S: es esperable que salary y yrs.service muestren curvatura en las colas (cola derecha más pesada), reflejando la asimetría observada en los histogramas.

5.3 Decisión sobre el tipo de correlación a utilizar

Regla de decisión adoptada en este análisis: si todas las variables involucradas en un par (la explicativa y salary) son normales según la prueba K-S, se utiliza el coeficiente de correlación de Pearson (mide asociación lineal). Si al menos una de las dos variables del par no es normal, se utiliza el coeficiente de correlación de Spearman (mide asociación monótona basada en rangos, y no requiere normalidad).

# Función de ayuda para decidir el tipo de correlación según el resultado de K-S
normal_phd     <- ks_phd$p.value     >= 0.05
normal_service <- ks_service$p.value >= 0.05
normal_salary  <- ks_salary$p.value  >= 0.05

cat("¿yrs.since.phd es normal? ", normal_phd, "\n")

## ¿yrs.since.phd es normal?  FALSE

cat("¿yrs.service es normal?   ", normal_service, "\n")

## ¿yrs.service es normal?    FALSE

cat("¿salary es normal?        ", normal_salary, "\n")

## ¿salary es normal?         FALSE

6 Análisis de correlación

6.1 Correlación entre variables numéricas y el salario

# Correlación de Spearman (rangos) - usada si alguna variable no es normal
cor_spearman_phd     <- cor.test(datos$yrs.since.phd, datos$salary, method = "spearman")
cor_spearman_service <- cor.test(datos$yrs.service, datos$salary, method = "spearman")

# Correlación de Pearson (lineal) - se calcula también como referencia
cor_pearson_phd     <- cor.test(datos$yrs.since.phd, datos$salary, method = "pearson")
cor_pearson_service <- cor.test(datos$yrs.service, datos$salary, method = "pearson")

cor_spearman_phd

## 
##  Spearman's rank correlation rho
## 
## data:  datos$yrs.since.phd and datos$salary
## S = 5468989, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4755676

cor_spearman_service

## 
##  Spearman's rank correlation rho
## 
## data:  datos$yrs.service and datos$salary
## S = 5996863, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4249486

tabla_cor <- data.frame(
  Variable = c("yrs.since.phd", "yrs.service"),
  Spearman_rho = c(cor_spearman_phd$estimate, cor_spearman_service$estimate),
  Spearman_p = c(cor_spearman_phd$p.value, cor_spearman_service$p.value),
  Pearson_r = c(cor_pearson_phd$estimate, cor_pearson_service$estimate),
  Pearson_p = c(cor_pearson_phd$p.value, cor_pearson_service$p.value)
) %>%
  mutate(
    Metodo_usado = "Spearman (por no cumplir normalidad)",
    Significativa_5pct = ifelse(Spearman_p < 0.05, "Sí", "No")
  )

kable(tabla_cor, digits = 4,
      caption = "Correlación entre variables numéricas explicativas y el salario") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Correlación entre variables numéricas explicativas y el salario
Variable	Spearman_rho	Spearman_p	Pearson_r	Pearson_p	Metodo_usado	Significativa_5pct
yrs.since.phd	0.4756	0	0.4192	0	Spearman (por no cumplir normalidad)	Sí
yrs.service	0.4249	0	0.3347	0	Spearman (por no cumplir normalidad)	Sí

6.1.1 Hipótesis de la prueba de correlación

\[ H_0: \rho = 0 \quad (\text{no existe correlación monótona entre la variable y el salario}) \] \[ H_1: \rho \neq 0 \quad (\text{existe correlación monótona entre la variable y el salario}) \]

Regla de decisión: se rechaza \(H_0\) si el valor p asociado es menor a 0.05.

Interpretación esperada: tanto yrs.since.phd como yrs.service deberían mostrar una correlación de Spearman positiva, moderada a alta, y estadísticamente significativa (valor p < 0.05) con salary. Esto confirma cuantitativamente lo observado en los diagramas de dispersión: a mayor antigüedad académica (medida de cualquiera de las dos formas), mayor es el salario. Como ambas variables resultan significativas, ambas se incluyen como candidatas para el modelo de regresión, conforme al criterio solicitado de incluir únicamente variables significativas.

6.2 Matriz de correlación visual

matriz_cor <- cor(datos %>% dplyr::select(yrs.since.phd, yrs.service, salary), method = "spearman")

kable(matriz_cor, digits = 3, caption = "Matriz de correlación de Spearman") %>%
  kable_styling(full_width = FALSE)

Matriz de correlación de Spearman
	yrs.since.phd	yrs.service	salary
yrs.since.phd	1.000	0.906	0.476
yrs.service	0.906	1.000	0.425
salary	0.476	0.425	1.000

# Visualización con corrplot si está disponible, si no, con un heatmap simple en ggplot
if (requireNamespace("corrplot", quietly = TRUE)) {
  corrplot::corrplot(matriz_cor, method = "color", type = "upper",
                      addCoef.col = "black", tl.col = "black", number.cex = 0.9)
} else {
  matriz_larga <- as.data.frame(as.table(matriz_cor))
  ggplot(matriz_larga, aes(Var1, Var2, fill = Freq)) +
    geom_tile() +
    geom_text(aes(label = round(Freq, 2)), color = "white") +
    scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B", midpoint = 0) +
    theme_minimal() +
    labs(title = "Matriz de correlación (Spearman)", x = "", y = "", fill = "rho")
}

Nota importante sobre yrs.since.phd y yrs.service: se observa (y es esperable) una correlación muy alta entre estas dos variables, ya que ambas miden, desde ángulos distintos, la antigüedad profesional de un docente. Esta alta correlación entre dos variables explicativas es justamente lo que en regresión se denomina multicolinealidad, un problema que se evaluará formalmente más adelante con el VIF y el GVIF, y que debe tenerse presente al construir los modelos.

6.3 Panel de correlaciones (Chart Correlation)

datos_numericos <- datos %>% dplyr::select(yrs.since.phd, yrs.service, salary)

chart.Correlation(datos_numericos, histogram = TRUE, method = "spearman", pch = 19,
                   col = "#AA96DA")

6.3.1 Interpretación del Chart Correlation y significancia de las variables

Lectura del panel:

salary y yrs.since.phd: se espera un coeficiente de correlación alto y positivo, acompañado de varios asteriscos (típicamente ***), lo que indica que la asociación entre los años desde el doctorado y el salario es estadísticamente significativa y muy poco probable que se deba al azar. Esto coincide con el resultado obtenido previamente mediante cor.test() con el método de Spearman.
salary y yrs.service: de manera similar, se espera una correlación positiva, significativa, aunque normalmente algo más baja que la anterior, ya que los años de servicio no siempre coinciden exactamente con los años desde el doctorado (un profesor puede haber trabajado en otra institución antes de llegar a la actual).
yrs.since.phd y yrs.service: esta es la correlación más alta del panel, y confirma visualmente (en la nube de puntos casi alineada sobre una recta creciente) el problema de multicolinealidad mencionado anteriormente. Dos variables tan correlacionadas entre sí aportan información muy similar al modelo, por lo que incluir ambas simultáneamente puede inflar la varianza de los coeficientes estimados, tal como se verificará formalmente en la sección de multicolinealidad (VIF y GVIF).

En síntesis, el panel de correlación confirma que ambas variables numéricas (yrs.since.phd y yrs.service) están significativamente asociadas con el salario, por lo que ambas son, en principio, candidatas válidas para el modelo de regresión, en concordancia con el criterio de incluir únicamente variables explicativas significativas. Sin embargo, su alta correlación mutua exige prestar especial atención a la multicolinealidad al construir los modelos múltiples.

7 Construcción del modelo de regresión lineal múltiple

7.1 Selección de variables para el modelo

A partir del análisis exploratorio y del análisis de correlación realizado, se concluyó que:

yrs.since.phd y yrs.service son variables numéricas significativamente correlacionadas con salary (correlación de Spearman, valor p < 0.05), por lo que ambas son candidatas válidas.
Las variables categóricas rank y discipline mostraron, en los boxplots, diferencias visibles en la mediana salarial entre sus grupos, por lo que también se consideran candidatas relevantes.
Existe alta colinealidad entre yrs.since.phd y yrs.service, lo cual deberá vigilarse en los modelos que incluyan ambas variables simultáneamente.

7.2 Modelo 1: salario en función de los años desde el doctorado y el rango académico

modelo1 <- lm(salary ~ yrs.since.phd + rank, data = datos)
summary(modelo1)

## 
## Call:
## lm(formula = salary ~ yrs.since.phd + rank, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -67148 -16683  -1323  11835 105552 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   81186.23    2964.17  27.389  < 2e-16 ***
## yrs.since.phd   -80.37     129.47  -0.621  0.53510    
## rankAssocProf 13932.18    4345.76   3.206  0.00146 ** 
## rankProf      47860.42    4412.63  10.846  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23650 on 393 degrees of freedom
## Multiple R-squared:  0.3948, Adjusted R-squared:  0.3902 
## F-statistic: 85.47 on 3 and 393 DF,  p-value: < 2.2e-16

7.2.1 Ecuación matemática del Modelo 1

Como rank es una variable categórica con 3 niveles (AsstProf, AssocProf, Prof), R la convierte automáticamente en 2 variables indicadoras (dummy), usando AsstProf como categoría de referencia. La ecuación general del modelo es:

\[ \widehat{salary}_i = \beta_0 + \beta_1 (\text{yrs.since.phd}_i) + \beta_2 (\text{rank\_AssocProf}_i) + \beta_3 (\text{rank\_Prof}_i) + \varepsilon_i \]

donde:

\(\text{rank\_AssocProf}_i = 1\) si el profesor \(i\) es Profesor Asociado, y 0 en caso contrario.
\(\text{rank\_Prof}_i = 1\) si el profesor \(i\) es Profesor Titular, y 0 en caso contrario.
Cuando ambas variables indicadoras son 0, el profesor pertenece a la categoría de referencia (AsstProf, Profesor Asistente).

7.2.2 Interpretación de los coeficientes del Modelo 1

\(\beta_0\) (intercepto): representa el salario promedio estimado para un profesor asistente (categoría de referencia) con 0 años desde el doctorado. Es un valor de referencia matemática, no siempre tiene sentido práctico directo si ningún profesor asistente tiene 0 años desde el doctorado en la muestra.
\(\beta_1\): indica el incremento promedio esperado en el salario por cada año adicional transcurrido desde el doctorado, manteniendo fijo el rango académico. Se espera un valor positivo, reflejando que la antigüedad académica se asocia con mayores ingresos.
\(\beta_2\): indica la diferencia promedio en salario entre un Profesor Asociado y un Profesor Asistente, controlando por los años desde el doctorado. Se espera un valor positivo.
\(\beta_3\): indica la diferencia promedio en salario entre un Profesor Titular y un Profesor Asistente, controlando por los años desde el doctorado. Se espera el valor positivo más alto de los tres coeficientes de rango, dado que Prof es el rango más alto.

7.3 Modelo 2: salario en función de los años de servicio y la disciplina

modelo2 <- lm(salary ~ yrs.service + discipline, data = datos)
summary(modelo2)

## 
## Call:
## lm(formula = salary ~ yrs.service + discipline, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -77537 -19699  -5135  15631 106625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  91335.8     3005.4  30.391  < 2e-16 ***
## yrs.service    862.8      109.2   7.904 2.73e-14 ***
## disciplineB  13184.0     2846.8   4.631 4.95e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27870 on 394 degrees of freedom
## Multiple R-squared:  0.1579, Adjusted R-squared:  0.1536 
## F-statistic: 36.94 on 2 and 394 DF,  p-value: 1.983e-15

7.3.1 Ecuación matemática del Modelo 2

Aquí discipline tiene 2 niveles (A y B), por lo que se genera una única variable indicadora, usando A como categoría de referencia:

\[ \widehat{salary}_i = \beta_0 + \beta_1 (\text{yrs.service}_i) + \beta_2 (\text{discipline\_B}_i) + \varepsilon_i \]

donde \(\text{discipline\_B}_i = 1\) si el profesor \(i\) pertenece a la disciplina B (aplicada), y 0 si pertenece a la disciplina A (teórica).

7.3.2 Interpretación de los coeficientes del Modelo 2

\(\beta_0\): salario promedio estimado para un profesor de disciplina A con 0 años de servicio.
\(\beta_1\): incremento promedio esperado en el salario por cada año adicional de servicio, manteniendo fija la disciplina.
\(\beta_2\): diferencia promedio en salario entre la disciplina B y la disciplina A, controlando por los años de servicio. Un valor positivo indicaría que, en promedio, los profesores de áreas aplicadas (B) ganan más que los de áreas teóricas (A), para el mismo número de años de servicio.

7.4 Modelo 3 (modelo completo propuesto): salario en función de ambas variables numéricas y ambas variables categóricas principales

modelo3 <- lm(salary ~ yrs.since.phd + yrs.service + rank + discipline, data = datos)
summary(modelo3)

## 
## Call:
## lm(formula = salary ~ yrs.since.phd + yrs.service + rank + discipline, 
##     data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65244 -13498  -1455   9638  99682 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    69869.0     3332.1  20.968  < 2e-16 ***
## yrs.since.phd    534.6      241.2   2.217  0.02720 *  
## yrs.service     -476.7      211.8  -2.250  0.02497 *  
## rankAssocProf  12831.5     4147.7   3.094  0.00212 ** 
## rankProf       45287.7     4236.7  10.689  < 2e-16 ***
## disciplineB    14505.2     2343.4   6.190 1.52e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22550 on 391 degrees of freedom
## Multiple R-squared:  0.4525, Adjusted R-squared:  0.4455 
## F-statistic: 64.64 on 5 and 391 DF,  p-value: < 2.2e-16

7.4.1 Ecuación matemática del Modelo 3

\[ \begin{aligned} \widehat{salary}_i = \; & \beta_0 + \beta_1 (\text{yrs.since.phd}_i) + \beta_2 (\text{yrs.service}_i) \\ & + \beta_3 (\text{rank\_AssocProf}_i) + \beta_4 (\text{rank\_Prof}_i) + \beta_5 (\text{discipline\_B}_i) + \varepsilon_i \end{aligned} \]

7.4.2 Interpretación de los coeficientes del Modelo 3

\(\beta_1\): efecto marginal de los años desde el doctorado sobre el salario, controlando simultáneamente por años de servicio, rango y disciplina.
\(\beta_2\): efecto marginal de los años de servicio sobre el salario, controlando por las demás variables. Es interesante notar que, al incluir ambas variables de antigüedad juntas, es posible que alguno de los dos coeficientes pierda significancia estadística o incluso cambie de signo respecto a los modelos simples, precisamente como consecuencia de la multicolinealidad detectada previamente entre ambas variables.
\(\beta_3, \beta_4\): diferencias salariales promedio de cada rango académico respecto a la categoría de referencia (AsstProf), controlando por antigüedad y disciplina.
\(\beta_5\): diferencia salarial promedio entre disciplinas, controlando por antigüedad y rango.

7.5 Comparación de los modelos propuestos

comparacion <- data.frame(
  Modelo = c("Modelo 1: yrs.since.phd + rank",
             "Modelo 2: yrs.service + discipline",
             "Modelo 3: yrs.since.phd + yrs.service + rank + discipline"),
  R2_ajustado = c(summary(modelo1)$adj.r.squared,
                  summary(modelo2)$adj.r.squared,
                  summary(modelo3)$adj.r.squared),
  AIC = c(AIC(modelo1), AIC(modelo2), AIC(modelo3)),
  BIC = c(BIC(modelo1), BIC(modelo2), BIC(modelo3))
)

kable(comparacion, digits = 4,
      caption = "Comparación de criterios de ajuste entre los tres modelos propuestos") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Comparación de criterios de ajuste entre los tres modelos propuestos
Modelo	R2_ajustado	AIC	BIC
Modelo 1: yrs.since.phd + rank	0.3902	9129.156	9149.076
Modelo 2: yrs.service + discipline	0.1536	9258.333	9274.269
Modelo 3: yrs.since.phd + yrs.service + rank + discipline	0.4455	9093.388	9121.275

Interpretación: el R cuadrado ajustado penaliza la inclusión de variables que no aportan capacidad explicativa real, por lo que es más adecuado que el R cuadrado simple para comparar modelos con distinto número de variables. Los criterios AIC y BIC (Criterio de Información de Akaike y de Schwarz/Bayesiano, respectivamente) penalizan también la complejidad del modelo; a menor valor de AIC o BIC, mejor es el balance entre ajuste y parsimonia. Se espera que el Modelo 3, al incluir más información relevante, tenga el mayor R cuadrado ajustado y los menores valores de AIC y BIC, por lo que se selecciona como el modelo final sobre el cual se realizará todo el diagnóstico de residuos a continuación.

modelo_final <- modelo3

8 Diagnóstico de residuos del modelo final

El diagnóstico de residuos es la etapa que permite verificar si se cumplen los supuestos del modelo de regresión lineal: linealidad, normalidad de los residuos, homocedasticidad (varianza constante) e independencia. Si estos supuestos no se cumplen, las pruebas de hipótesis sobre los coeficientes (valores p, intervalos de confianza) pueden no ser confiables.

8.1 Extracción de los residuos

residuos <- residuals(modelo_final)
residuos_estandarizados <- rstandard(modelo_final)
residuos_estudentizados <- rstudent(modelo_final)
valores_ajustados <- fitted(modelo_final)

8.2 Prueba de normalidad de los residuos (Kolmogorov-Smirnov)

Tal como se hizo con las variables originales, y dado que el tamaño de la muestra es n = 397 (n > 50), se utiliza nuevamente la prueba de Kolmogorov-Smirnov para evaluar si los residuos del modelo siguen una distribución normal.

8.2.1 Planteamiento de hipótesis

\[ H_0: \text{Los residuos del modelo provienen de una distribución normal} \] \[ H_1: \text{Los residuos del modelo NO provienen de una distribución normal} \]

ks_residuos <- ks.test(scale(residuos), "pnorm")
ks_residuos

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  scale(residuos)
## D = 0.090618, p-value = 0.002947
## alternative hypothesis: two-sided

Regla de decisión: si el valor p obtenido es menor a 0.05, se rechaza \(H_0\) y se concluye que los residuos no son normales, lo cual sugeriría que el modelo aún no captura toda la estructura de los datos (posible falta de linealidad, variables omitidas, o presencia de valores atípicos influyentes que distorsionan el ajuste).

8.3 Gráfico Q-Q de los residuos

qqnorm(residuos_estandarizados, main = "Q-Q Plot de los residuos estandarizados",
       col = "#AA96DA", pch = 19)
qqline(residuos_estandarizados, col = "#D7191C", lwd = 2)

Interpretación: si los puntos se ajustan razonablemente a la línea diagonal roja, especialmente en la parte central, se considera que el supuesto de normalidad es aceptable. Desviaciones en los extremos (colas) suelen indicar la presencia de valores atípicos, los cuales se analizarán formalmente más adelante con la prueba de Bonferroni.

8.4 Análisis gráfico global del modelo: plot del summary (gráficos de diagnóstico de `lm`)

R permite generar automáticamente cuatro gráficos de diagnóstico estándar a partir de un objeto lm, los cuales en conjunto permiten validar visualmente los residuos del modelo.

par(mfrow = c(2, 2))
plot(modelo_final, col = "#AA96DA", pch = 19)

par(mfrow = c(1, 1))

Interpretación de cada panel:

Residuals vs Fitted (Residuos vs. valores ajustados): permite evaluar la linealidad y la homocedasticidad. Si los puntos se distribuyen aleatoriamente alrededor de la línea horizontal en cero, sin un patrón claro (forma de embudo, curva, etc.), los supuestos de linealidad y varianza constante son razonables.
Q-Q Residuals: el mismo gráfico Q-Q visto anteriormente, generado automáticamente por R, para verificar normalidad.
Scale-Location: muestra la raíz cuadrada de los residuos estandarizados en valor absoluto contra los valores ajustados. Una línea de tendencia aproximadamente horizontal apoya el supuesto de homocedasticidad; una tendencia creciente sugeriría heterocedasticidad (la varianza del error aumenta con el nivel del salario ajustado).
Residuals vs Leverage: permite identificar observaciones con alto apalancamiento (leverage) y su posible influencia sobre el modelo, mostrando además las curvas de la Distancia de Cook, que se retomará formalmente más adelante.

8.5 Prueba de homocedasticidad (varianza constante de los residuos)

Para complementar el análisis visual del gráfico Scale-Location, se aplica formalmente la prueba de Breusch-Pagan, ampliamente utilizada para contrastar la homocedasticidad en modelos de regresión lineal.

8.5.1 Planteamiento de hipótesis

\[ H_0: \text{La varianza de los residuos es constante (homocedasticidad)} \] \[ H_1: \text{La varianza de los residuos NO es constante (heterocedasticidad)} \]

bp_test <- bptest(modelo_final)
bp_test

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_final
## BP = 63.172, df = 5, p-value = 2.681e-12

Regla de decisión: si el valor p es menor a 0.05, se rechaza \(H_0\) y se concluye que existe heterocedasticidad, es decir, que la dispersión de los errores no es constante a lo largo de los valores ajustados del salario (por ejemplo, podría haber mayor variabilidad en los salarios altos que en los bajos). Si esto ocurriera, se recomendaría considerar transformaciones de la variable respuesta (como el logaritmo del salario) o el uso de errores estándar robustos.

8.6 Prueba de independencia de los residuos (Durbin-Watson)

La independencia de los residuos es relevante principalmente cuando existe un orden natural en las observaciones (por ejemplo, series de tiempo). Aunque en este conjunto de datos las observaciones no tienen un orden temporal explícito, se incluye la prueba de Durbin-Watson por ser la prueba estándar de independencia (autocorrelación) de los residuos, tal como fue solicitado.

8.6.1 Planteamiento de hipótesis

\[ H_0: \text{No existe autocorrelación entre los residuos (los residuos son independientes)} \] \[ H_1: \text{Existe autocorrelación entre los residuos} \]

dw_test <- durbinWatsonTest(modelo_final)
dw_test

##  lag Autocorrelation D-W Statistic p-value
##    1      0.04953982      1.900167   0.298
##  Alternative hypothesis: rho != 0

Regla de decisión e interpretación: el estadístico de Durbin-Watson toma valores entre 0 y 4. Un valor cercano a 2 indica ausencia de autocorrelación (no se rechaza \(H_0\)); valores cercanos a 0 sugieren autocorrelación positiva, y valores cercanos a 4 sugieren autocorrelación negativa. Adicionalmente, si el valor p reportado es mayor a 0.05, se concluye que no hay evidencia suficiente para rechazar la independencia de los residuos.

9 Detección de valores atípicos (outliers): prueba de Bonferroni

Un outlier (valor atípico) es una observación cuyo residuo es inusualmente grande, es decir, cuyo valor real del salario se aleja mucho de lo que el modelo predice para esa observación. La función outlierTest() del paquete car identifica las observaciones con el residuo estudentizado más extremo y aplica una corrección de Bonferroni sobre el valor p, debido a que se están realizando múltiples comparaciones simultáneas (una por cada observación de la muestra), lo cual incrementa la probabilidad de cometer un error de tipo I si no se corrige.

9.0.1 Planteamiento de hipótesis

\[ H_0: \text{La observación NO es un valor atípico (su residuo proviene de la misma distribución que el resto)} \] \[ H_1: \text{La observación SÍ es un valor atípico} \]

outlier_test <- outlierTest(modelo_final)
outlier_test

##     rstudent unadjusted p-value Bonferroni p
## 44  4.559584         6.8711e-06    0.0027278
## 365 3.999664         7.5895e-05    0.0301300

Regla de decisión: se considera que una observación es un outlier estadísticamente significativo si su valor p ajustado por Bonferroni (columna Bonferroni p) es menor a 0.05. Si el resultado no muestra observaciones significativas, se concluye que, aunque puedan existir residuos grandes en términos absolutos, ninguno es lo suficientemente extremo como para considerarse un valor atípico estadísticamente comprobado, una vez corregido por el número de comparaciones realizadas.

plot(residuos_estudentizados, pch = 19, col = "#AA96DA",
     main = "Residuos estudentizados por observación",
     xlab = "Índice de observación", ylab = "Residuo estudentizado")
abline(h = c(-2, 2), col = "#D7191C", lty = 2, lwd = 1.5)
abline(h = c(-3, 3), col = "#6A4C93", lty = 2, lwd = 1.5)
legend("topright", legend = c("Límite ±2", "Límite ±3"),
       col = c("#D7191C", "#6A4C93"), lty = 2, cex = 0.8)

Interpretación: como regla práctica complementaria, los residuos estudentizados con magnitud mayor a 2 (en valor absoluto) se consideran observaciones inusuales, y aquellos mayores a 3 se consideran claramente atípicos. Las líneas de referencia en el gráfico permiten identificar visualmente cuántas observaciones superan estos umbrales.

10 Puntos de apalancamiento: matriz sombrero (Hat Matrix)

Un punto de apalancamiento (leverage point) es una observación cuyos valores en las variables explicativas (X) son inusuales o extremos, independientemente de si su valor de Y (salario) está bien o mal predicho. Estos puntos tienen el potencial de “halar” la recta de regresión hacia sí mismos.

El apalancamiento de cada observación se obtiene de la diagonal de la matriz sombrero (Hat Matrix), definida como:

\[ H = X (X^T X)^{-1} X^T \]

donde \(X\) es la matriz de diseño del modelo (incluye el intercepto y todas las variables explicativas, ya codificadas). Los valores de la diagonal de \(H\), denotados \(h_{ii}\), cumplen \(0 \leq h_{ii} \leq 1\), y la suma de todos los \(h_{ii}\) es igual al número de parámetros estimados en el modelo (\(p\), incluyendo el intercepto).

valores_hat <- hatvalues(modelo_final)

# Número de parámetros (p) y número de observaciones (n)
p <- length(coef(modelo_final))
n <- nrow(datos)

# Regla práctica de corte para identificar alto apalancamiento: 2p/n o 3p/n
corte_leverage_1 <- 2 * p / n
corte_leverage_2 <- 3 * p / n

cat("Número de parámetros del modelo (p):", p, "\n")

## Número de parámetros del modelo (p): 6

cat("Número de observaciones (n):", n, "\n")

## Número de observaciones (n): 397

cat("Punto de corte 2p/n:", round(corte_leverage_1, 4), "\n")

## Punto de corte 2p/n: 0.0302

cat("Punto de corte 3p/n:", round(corte_leverage_2, 4), "\n")

## Punto de corte 3p/n: 0.0453

puntos_alto_leverage <- which(valores_hat > corte_leverage_2)
cat("Número de observaciones con alto apalancamiento (> 3p/n):", length(puntos_alto_leverage), "\n")

## Número de observaciones con alto apalancamiento (> 3p/n): 4

plot(valores_hat, pch = 19, col = "#A0CED9",
     main = "Valores de apalancamiento (hat values) por observación",
     xlab = "Índice de observación", ylab = "Valor hat (h_ii)")
abline(h = corte_leverage_1, col = "#D7191C", lty = 2, lwd = 1.5)
abline(h = corte_leverage_2, col = "#6A4C93", lty = 2, lwd = 1.5)
legend("topright", legend = c("Corte 2p/n", "Corte 3p/n"),
       col = c("#D7191C", "#6A4C93"), lty = 2, cex = 0.8)

Interpretación: las observaciones cuyo valor \(h_{ii}\) supera el punto de corte (usualmente \(2p/n\) o, de forma más conservadora, \(3p/n\)) se consideran puntos de alto apalancamiento, es decir, tienen combinaciones de valores en las variables explicativas (por ejemplo, una antigüedad extremadamente alta o baja) poco frecuentes dentro de la muestra. Es importante recalcar que un punto de apalancamiento no necesariamente es un outlier ni un punto influyente: solo indica que su posición en el espacio de las variables X es atípica; su efecto real sobre el modelo se evalúa con las medidas de influencia que se presentan a continuación.

11 Puntos influyentes: Distancia de Cook, DFFITS y DFBETAS

Un punto influyente es aquella observación que, al ser eliminada de la muestra, cambia de forma notable las estimaciones del modelo (los coeficientes, los valores ajustados, o ambos). Un punto influyente combina, en cierta medida, ser un outlier (residuo grande) y tener cierto apalancamiento.

11.1 Distancia de Cook

La Distancia de Cook mide cuánto cambiarían todos los valores ajustados del modelo si se eliminara la observación \(i\)-ésima. Se calcula como:

\[ D_i = \frac{(\hat{\beta} - \hat{\beta}_{(i)})^T (X^T X) (\hat{\beta} - \hat{\beta}_{(i)})}{p \cdot S^2} \]

donde \(\hat{\beta}_{(i)}\) son los coeficientes estimados al excluir la observación \(i\), \(p\) es el número de parámetros del modelo y \(S^2\) es la varianza estimada de los residuos.

distancias_cook <- cooks.distance(modelo_final)

# Punto de corte habitual: 4/n (regla práctica) y también 1 (regla clásica más conservadora)
corte_cook <- 4 / n

plot(distancias_cook, pch = 19, col = "#FBC4AB",
     main = "Distancia de Cook por observación",
     xlab = "Índice de observación", ylab = "Distancia de Cook")
abline(h = corte_cook, col = "#D7191C", lty = 2, lwd = 1.5)
legend("topright", legend = paste0("Corte 4/n = ", round(corte_cook, 4)),
       col = "#D7191C", lty = 2, cex = 0.8)

puntos_influyentes_cook <- which(distancias_cook > corte_cook)
cat("Número de observaciones influyentes según Distancia de Cook (> 4/n):",
    length(puntos_influyentes_cook), "\n")

## Número de observaciones influyentes según Distancia de Cook (> 4/n): 18

Interpretación: las observaciones cuya Distancia de Cook supera el punto de corte (\(4/n\), una regla práctica comúnmente usada; algunos autores también usan el valor 1 como corte más conservador) se consideran potencialmente influyentes y merecen revisión adicional (verificar si son errores de registro, casos genuinamente atípicos, o simplemente observaciones legítimas en los extremos de la distribución).

11.2 DFFITS

El estadístico DFFITS mide el cambio en el valor ajustado de la propia observación \(i\) (no de todo el modelo, como la Distancia de Cook) al excluirla del ajuste, estandarizado por su error estándar:

\[ DFFITS_i = \frac{\hat{y}_i - \hat{y}_{i(i)}}{s_{(i)} \sqrt{h_{ii}}} \]

valores_dffits <- dffits(modelo_final)

# Punto de corte habitual: 2 * sqrt(p/n)
corte_dffits <- 2 * sqrt(p / n)

plot(valores_dffits, pch = 19, col = "#B5EAD7",
     main = "DFFITS por observación",
     xlab = "Índice de observación", ylab = "DFFITS")
abline(h = c(-corte_dffits, corte_dffits), col = "#D7191C", lty = 2, lwd = 1.5)
legend("topright", legend = paste0("Corte ±2*sqrt(p/n) = ±", round(corte_dffits, 4)),
       col = "#D7191C", lty = 2, cex = 0.8)

puntos_influyentes_dffits <- which(abs(valores_dffits) > corte_dffits)
cat("Número de observaciones influyentes según DFFITS:",
    length(puntos_influyentes_dffits), "\n")

## Número de observaciones influyentes según DFFITS: 19

Interpretación: un valor absoluto de DFFITS mayor al punto de corte \(2\sqrt{p/n}\) indica que esa observación, individualmente, tiene una influencia considerable sobre su propio valor predicho. Es un complemento útil a la Distancia de Cook porque se enfoca en el efecto sobre la predicción individual, en lugar del efecto sobre el conjunto completo de coeficientes.

11.3 DFBETAS

El estadístico DFBETAS mide el cambio, en unidades de error estándar, en cada coeficiente individual del modelo al excluir la observación \(i\). Se calcula uno por cada combinación de observación y coeficiente, lo que permite identificar no solo qué observaciones son influyentes, sino sobre cuál coeficiente específico ejercen mayor influencia.

valores_dfbetas <- dfbetas(modelo_final)

# Punto de corte habitual: 2/sqrt(n)
corte_dfbetas <- 2 / sqrt(n)

# Visualización para el coeficiente de yrs.since.phd como ejemplo
plot(valores_dfbetas[, "yrs.since.phd"], pch = 19, col = "#FFDAC1",
     main = "DFBETAS para el coeficiente de yrs.since.phd",
     xlab = "Índice de observación", ylab = "DFBETAS (yrs.since.phd)")
abline(h = c(-corte_dfbetas, corte_dfbetas), col = "#D7191C", lty = 2, lwd = 1.5)
legend("topright", legend = paste0("Corte ±2/sqrt(n) = ±", round(corte_dfbetas, 4)),
       col = "#D7191C", lty = 2, cex = 0.8)

# Conteo de observaciones influyentes por cada coeficiente del modelo
conteo_dfbetas <- colSums(abs(valores_dfbetas) > corte_dfbetas)

kable(data.frame(Coeficiente = names(conteo_dfbetas),
                  Num_observaciones_influyentes = conteo_dfbetas),
      caption = "Número de observaciones influyentes por coeficiente, según DFBETAS") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Número de observaciones influyentes por coeficiente, según DFBETAS
	Coeficiente	Num_observaciones_influyentes
(Intercept)	(Intercept)	4
yrs.since.phd	yrs.since.phd	26
yrs.service	yrs.service	29
rankAssocProf	rankAssocProf	4
rankProf	rankProf	10
disciplineB	disciplineB	21

Interpretación: un valor absoluto de DFBETAS mayor al corte \(2/\sqrt{n}\) para un coeficiente específico indica que esa observación, de ser eliminada, cambiaría notablemente la estimación de ese coeficiente en particular. Esto es especialmente útil para detectar, por ejemplo, si un solo profesor con un salario o una antigüedad muy particular está “empujando” artificialmente el coeficiente asociado a yrs.since.phd o a algún nivel de rank.

11.4 Resumen gráfico conjunto de influencia

influencePlot(modelo_final, col = "#AA96DA",
              main = "Gráfico de influencia: leverage, residuos estudentizados y Distancia de Cook")

Interpretación: este gráfico, generado con influencePlot() del paquete car, combina en un solo lugar los tres conceptos trabajados: en el eje horizontal se ubica el apalancamiento (hat values), en el eje vertical los residuos estudentizados, y el tamaño del círculo de cada punto es proporcional a su Distancia de Cook. Las observaciones que se destacan por estar muy a la derecha (alto apalancamiento), muy arriba o abajo (residuo grande), o con un círculo grande (alta Distancia de Cook) son las que merecen mayor atención como posibles puntos problemáticos para el modelo.

12 Análisis de multicolinealidad: VIF y GVIF

La multicolinealidad ocurre cuando dos o más variables explicativas del modelo están altamente correlacionadas entre sí, lo que dificulta separar el efecto individual de cada una sobre la variable respuesta, e infla la varianza (y por tanto el error estándar) de los coeficientes estimados, restando precisión y estabilidad al modelo.

12.1 Factor de Inflación de la Varianza (VIF)

Para variables explicativas numéricas o categóricas de solo 2 niveles, el VIF de la variable \(X_k\) se calcula como:

\[ VIF_k = \frac{1}{1 - R_k^2} \]

donde \(R_k^2\) es el coeficiente de determinación obtenido al regresar la variable \(X_k\) contra todas las demás variables explicativas del modelo.

12.2 Factor de Inflación de Varianza Generalizado (GVIF)

Cuando una variable categórica tiene más de 2 niveles (como rank, con 3 niveles), se representa mediante varias variables indicadoras (dummy) simultáneamente, por lo que el VIF tradicional (pensado para una sola columna) no aplica directamente. En su lugar, se utiliza el GVIF, propuesto por Fox y Monette (1992), que generaliza el concepto de inflación de varianza para conjuntos de variables relacionadas (como todas las dummies de un mismo factor categórico). Para poder comparar de forma justa variables con distinto número de grados de libertad (df), se suele reportar también:

\[ GVIF^{1/(2 \cdot df)} \]

Este valor ajustado es comparable, en magnitud, al \(\sqrt{VIF}\) de una variable numérica tradicional, lo que facilita la interpretación conjunta.

vif_resultado <- vif(modelo_final)
vif_resultado

##                   GVIF Df GVIF^(1/(2*Df))
## yrs.since.phd 7.518920  1        2.742065
## yrs.service   5.908984  1        2.430840
## rank          2.003562  2        1.189736
## discipline    1.063139  1        1.031086

# La función vif() de car devuelve GVIF automáticamente si hay factores con más de 2 niveles
if (is.matrix(vif_resultado)) {
  tabla_vif <- as.data.frame(vif_resultado)
  tabla_vif$Variable <- rownames(tabla_vif)
  rownames(tabla_vif) <- NULL
  tabla_vif <- tabla_vif[, c("Variable", "GVIF", "Df", "GVIF^(1/(2*Df))")]
} else {
  tabla_vif <- data.frame(Variable = names(vif_resultado), VIF = vif_resultado)
}

kable(tabla_vif, digits = 3,
      caption = "Factor de Inflación de Varianza (VIF) / Generalizado (GVIF) del modelo final") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Factor de Inflación de Varianza (VIF) / Generalizado (GVIF) del modelo final
Variable	GVIF	Df	GVIF^(1/(2*Df))
yrs.since.phd	7.519	1	2.742
yrs.service	5.909	1	2.431
rank	2.004	2	1.190
discipline	1.063	1	1.031

12.2.1 Criterio de interpretación

La regla práctica más utilizada en la literatura estadística para interpretar el VIF (o el GVIF ajustado, \(GVIF^{1/(2 \cdot df)}\), elevado al cuadrado para comparar en la misma escala que el VIF) es la siguiente:

Valor de VIF (o GVIF ajustado al cuadrado)	Interpretación
VIF = 1	Ausencia total de correlación con las demás variables
1 < VIF < 5	Multicolinealidad moderada, generalmente aceptable
VIF ≥ 5 (algunos autores usan el corte en 10)	Multicolinealidad alta, preocupante para el modelo

Interpretación esperada para este modelo: dado que se identificó previamente una correlación muy alta entre yrs.since.phd y yrs.service, es razonable esperar que estas dos variables presenten los valores de VIF más elevados del modelo, posiblemente superando el umbral de 5. Si este fuera el caso, se recomendaría, como estrategia de mitigación, retirar una de las dos variables redundantes del modelo (por ejemplo, conservar solo yrs.since.phd, que captura de forma más directa la antigüedad académica desde el doctorado), o bien combinarlas en un solo índice, o aplicar técnicas de regresión penalizada (como Ridge) que no forman parte del alcance de este documento pero se mencionan como posible extensión del análisis. Las variables categóricas (rank, discipline), al no estar mecánicamente ligadas a las variables numéricas, deberían mostrar valores de GVIF ajustado cercanos a 1, indicando que no representan un problema de multicolinealidad relevante en este modelo.

13 Conclusiones generales

A partir de todo el análisis realizado sobre el conjunto de datos Salaries, se pueden extraer las siguientes conclusiones generales, las cuales deben ser confirmadas con los valores numéricos exactos obtenidos al ejecutar este documento en RStudio:

El salario de los profesores universitarios está asociado de forma significativa tanto con variables numéricas de antigüedad (yrs.since.phd, yrs.service) como con variables categóricas de jerarquía y área (rank, discipline).
Ninguna de las variables numéricas analizadas cumplió el supuesto de normalidad según la prueba de Kolmogorov-Smirnov, lo cual justificó el uso del coeficiente de correlación de Spearman en lugar de Pearson para evaluar su asociación con el salario.
El modelo de regresión lineal múltiple que combina ambas variables numéricas y ambas variables categóricas principales (Modelo 3) mostró el mejor desempeño relativo según el R cuadrado ajustado, el AIC y el BIC, por lo que se seleccionó como modelo final para el diagnóstico completo.
El diagnóstico de residuos (normalidad, homocedasticidad e independencia) permite evaluar la validez de las pruebas de hipótesis sobre los coeficientes del modelo; cualquier violación detectada (por ejemplo, heterocedasticidad) debe tenerse en cuenta al comunicar los resultados, y podría motivar el uso de transformaciones de la variable respuesta o de errores estándar robustos.
El análisis de valores atípicos, apalancamiento e influencia (outlierTest, hat values, Distancia de Cook, DFFITS y DFBETAS) permite identificar observaciones particulares que podrían estar afectando de manera desproporcionada las conclusiones del modelo, y que merecerían una revisión adicional antes de presentar el modelo como definitivo.
La alta correlación entre yrs.since.phd y yrs.service, confirmada cuantitativamente mediante el VIF y el GVIF, constituye la principal limitación del modelo completo y debe comunicarse con transparencia, junto con las posibles estrategias de mitigación discutidas.

Este documento integra, de manera secuencial y justificada estadísticamente, las etapas de exploración, prueba de hipótesis de normalidad, análisis de correlación, modelado, y diagnóstico de regresión lineal múltiple, ofreciendo una base sólida para el aprendizaje de estos conceptos en Estadística Aplicada.

Análisis Exploratorio y Modelo de Regresión Lineal Múltiple sobre el Conjunto de Datos Salaries

Maira Alejandra Vargas Castro

26 de junio de 2026