1 Objetivo y datos

Antes de analizar nada hay que saber qué representa cada fila, qué variables hay y cómo se relacionan las tablas. En este informe se trabaja con tres fuentes: iris (150 flores), dos extractos históricos ZA8630 (familias e hijos) y business_staff (1000 personas, corte transversal). Los datos se leen desde ../data, por lo que el .Rmd debe estar en la carpeta rmarkdown/.

library(data.table)
library(dplyr)
library(ggplot2)
library(lubridate)
library(modelsummary)
library(knitr)
iris_data      <- fread(file.path(data_path, "iris_data.txt"))
business_staff <- fread(file.path(data_path, "business_staff_cross_section_original_names.txt"))
families_raw   <- fread(file.path(data_path, "ZA8630_families_teaching_extract_original_names.txt"),
                        na.strings = c("", "NA"))
children_raw   <- fread(file.path(data_path, "ZA8630_children_teaching_extract_original_names.txt"),
                        na.strings = c("", "NA"))

Qué representa cada fila en cada tabla:

Tabla Unidad de observación (una fila = …) Filas
iris una flor 150
families (ZA8630) una familia (matrimonio) 34.708
children (ZA8630) un hijo 80.486
business_staff una persona (corte transversal) 1.000

2 Preparación y estructura de los datos

2.1 Decisiones de preparación de datos

Los archivos originales no se modifican: todos los cambios se hacen en R, sobre copias en memoria, de modo que los datos originales y los cálculos quedan separados (Broman y Woo, 2018).

# Nombres de trabajo (snake_case, en inglés)
setnames(business_staff,
         old = c("obs", "EDU", "inc", "zufr", "trip"),
         new = c("wave", "education", "income", "happiness", "trip_abroad"))

setnames(families_raw,
         old = c("doc", "dat4", "cort", "sozial", "grasen",
                 "idm", "gebm4", "todm4", "idf", "gebf4", "todf4"),
         new = c("family_id", "marriage_date", "marriage_place", "social_group",
                 "landholding_grasen", "husband_id", "husband_birth_date",
                 "husband_death_date", "wife_id", "wife_birth_date", "wife_death_date"))

setnames(children_raw,
         old = c("doc", "idk", "sex", "gebk4", "todk4", "gebortk", "todortk", "berufek"),
         new = c("family_id", "child_id", "gender", "child_birth_date",
                 "child_death_date", "child_birth_place", "child_death_place",
                 "child_occupation"))

# Recodificación de sexo: M/W -> male/female
children_raw[gender == "M", gender := "male"]
children_raw[gender == "W", gender := "female"]

# Fechas como Date
date_variables_families <- c("marriage_date", "husband_birth_date", "husband_death_date", "wife_birth_date", "wife_death_date")
families_raw[, (date_variables_families) := lapply(.SD, as.Date), .SDcols = date_variables_families]
children_raw[, child_birth_date := as.Date(child_birth_date)]
children_raw[, child_death_date := as.Date(child_death_date)]


# Factores en business_staff (low es la categoría de referencia)
business_staff[, education := factor(education, levels = c("low", "middle", "high"))]
business_staff[, gender    := factor(gender)]

Las decisiones principales y su justificación son:

  1. Nombres descriptivos y homogéneos. Los nombres técnicos (doc, gebk4, zufr, inc…) se sustituyen por nombres en inglés con un formato único (family_id, child_birth_date, happiness, income…). Broman y Woo (2018) recomiendan nombres consistentes, sin ambigüedad, y documentar la equivalencia en un diccionario de datos (ver Ejercicio 3).
  2. Fechas en formato Date (ISO 8601). Las fechas se leen como Date para poder restarlas y calcular edades; esto sigue la recomendación de Broman y Woo sobre fechas sin ambigüedad.
  3. Valores ausentes explícitos. Las cadenas vacías se leen como NA (na.strings = c("", "NA")), de modo que cada celda contiene un único valor y los ausentes se codifican de una sola manera.
  4. Recodificación de gender. M/W pasan a male/female, la misma convención que business_staff, lo que evita mezclar códigos en distintos archivos.
  5. education como factor ordenado por niveles. Los niveles low < middle < high hacen que las tablas y gráficos respeten el orden natural. Consecuencia: low es la categoría de referencia en las regresiones, y todos los coeficientes de education se leen como diferencias respecto a la educación baja.
  6. Una tabla por tipo de unidad observada. families (una fila por familia) y children (una fila por hijo) se mantienen separadas y solo se combinan para el análisis, según el principio de Wickham (2014).
  7. Variables derivadas aparte de los datos originales. razon, diferencia y las edades parentales se calculan en R, sin sobrescribir las columnas de partida.

2.2 Ejercicio 3. Archivos históricos ZA8630: correspondencia de nombres

families_raw contiene una fila por familia (matrimonio) y children_raw una fila por hijo. El script base ya renombró las variables; la equivalencia queda documentada así (un pequeño data dictionary, como recomiendan Broman y Woo, 2018).

dic_families <- data.frame(
  original = c("doc", "dat4", "cort", "sozial", "grasen", "idm", "gebm4",
               "todm4", "idf", "gebf4", "todf4"),
  trabajo  = c("family_id", "marriage_date", "marriage_place", "social_group",
               "landholding_grasen", "husband_id", "husband_birth_date",
               "husband_death_date", "wife_id", "wife_birth_date", "wife_death_date"),
  significado = c("Identificador de la familia", "Fecha del matrimonio",
                  "Lugar del matrimonio", "Grupo social",
                  "Tenencia de tierras (Grasen)", "Identificador del marido",
                  "Nacimiento del marido", "Fallecimiento del marido",
                  "Identificador de la esposa", "Nacimiento de la esposa",
                  "Fallecimiento de la esposa"))
kable(dic_families, caption = "Correspondencia de nombres: families")
Correspondencia de nombres: families
original trabajo significado
doc family_id Identificador de la familia
dat4 marriage_date Fecha del matrimonio
cort marriage_place Lugar del matrimonio
sozial social_group Grupo social
grasen landholding_grasen Tenencia de tierras (Grasen)
idm husband_id Identificador del marido
gebm4 husband_birth_date Nacimiento del marido
todm4 husband_death_date Fallecimiento del marido
idf wife_id Identificador de la esposa
gebf4 wife_birth_date Nacimiento de la esposa
todf4 wife_death_date Fallecimiento de la esposa 
dic_children <- data.frame(
  original = c("doc", "idk", "sex", "gebk4", "todk4", "gebortk", "todortk", "berufek"),
  trabajo  = c("family_id", "child_id", "gender", "child_birth_date",
               "child_death_date", "child_birth_place", "child_death_place",
               "child_occupation"),
  significado = c("Familia a la que pertenece el hijo", "Identificador del hijo",
                  "Sexo (M = male, W = female)", "Nacimiento del hijo",
                  "Fallecimiento del hijo", "Lugar de nacimiento",
                  "Lugar de fallecimiento", "Ocupación"))
kable(dic_children, caption = "Correspondencia de nombres: children")
Correspondencia de nombres: children
original trabajo significado
doc family_id Familia a la que pertenece el hijo
idk child_id Identificador del hijo
sex gender Sexo (M = male, W = female)
gebk4 child_birth_date Nacimiento del hijo
todk4 child_death_date Fallecimiento del hijo
gebortk child_birth_place Lugar de nacimiento
todortk child_death_place Lugar de fallecimiento
berufek child_occupation Ocupación

Además, sex se recodificó de M/W a male/female. La clave que une ambas tablas es family_id (antes doc en los dos archivos).

2.3 Ejercicio 4. Combinación de tablas y edades parentales

2.3.1 Por qué la relación es uno-a-muchos

anyDuplicated(families_raw$family_id)                       # 0: clave única en families
## [1] 0
all(children_raw$family_id %in% families_raw$family_id)     # todo hijo tiene familia
## [1] TRUE
hijos_por_familia <- children_raw[, .(n_children = .N), by = family_id]
sum(hijos_por_familia$n_children > 1)
## [1] 14756

family_id no se repite en families_raw, pero sí en children_raw: una familia puede tener varios hijos y cada hijo pertenece a una sola familia. Por eso la relación es uno-a-muchos. De hecho, 14.756 familias tienen más de un hijo.

2.3.2 Fusión y cálculo de edades

Se parte de children_raw con all.x = TRUE para conservar todos los hijos. Se descartan za_nr, version y doi (metadatos repetidos en ambas tablas) para evitar columnas duplicadas con sufijos .x/.y.

kid_fam <- merge(
  children_raw[, .(family_id, child_id, gender, child_birth_date)],
  families_raw[, .(family_id, marriage_date, husband_id, husband_birth_date,
                   wife_id, wife_birth_date)],
  by = "family_id", all.x = TRUE)

c(hijos_originales = nrow(children_raw), filas_tras_fusion = nrow(kid_fam))
##  hijos_originales filas_tras_fusion 
##             80486             80486

La fusión no pierde ni duplica hijos. Ahora se calculan las edades en años (diferencia en días entre fechas dividida por 365,25):

kid_fam[, `:=`(
  mother_age_at_birth = as.numeric(child_birth_date - wife_birth_date)    / 365.25,
  father_age_at_birth = as.numeric(child_birth_date - husband_birth_date) / 365.25)]

edad_media <- kid_fam[, .(
  n_madre_valida    = sum(!is.na(mother_age_at_birth)),
  edad_media_madres = mean(mother_age_at_birth, na.rm = TRUE),
  n_padre_valida    = sum(!is.na(father_age_at_birth)),
  edad_media_padres = mean(father_age_at_birth, na.rm = TRUE))]
kable(edad_media, digits = 2, caption = "Edad media al nacimiento del hijo")
Edad media al nacimiento del hijo
n_madre_valida edad_media_madres n_padre_valida edad_media_padres
52820 32.02 52348 35.72 
kid_fam[, .(min_madre = min(mother_age_at_birth, na.rm = TRUE),
            max_madre = max(mother_age_at_birth, na.rm = TRUE),
            min_padre = min(father_age_at_birth, na.rm = TRUE),
            max_padre = max(father_age_at_birth, na.rm = TRUE),
            edades_negativas = sum(mother_age_at_birth < 0 | father_age_at_birth < 0,
                                   na.rm = TRUE))]
##    min_madre max_madre min_padre max_padre edades_negativas
##        <num>     <num>     <num>     <num>            <int>
## 1: 0.2354552  55.77276  16.07392  78.39014                0

La edad media de las madres al nacimiento es de 32.02 años y la de los padres de 35.72, es decir, unos 3.7 años más. Las medias solo usan los casos con ambas fechas conocidas (52.820 en madres y 52.348 en padres, de 80.486 hijos), así que describen a las familias con fechas registradas, que pueden no ser representativas del conjunto.

2.3.3 Unidad de observación y tidy data

La unidad de observación de kid_fam es el hijo: cada fila es un hijo, cada columna una variable y cada celda un valor, de modo que cumple los criterios de tidy data de Wickham (2014). Los datos de la familia (fechas de los padres) se repiten para hermanos, pero esa redundancia no rompe el orden de la tabla; solo hay que tenerla en cuenta, por ejemplo, al contar familias con una tabla donde la unidad es el hijo.

2.3.4 Añadir las edades parentales a families

Esta operación no es directa: families tiene una fila por familia, mientras que las edades se calculan por hijo y pueden variar entre hermanos (los padres envejecen). Hay que decidir cómo resumir varias edades en un único valor por familia. Aquí se usa la media por familia, antes de unir; con otros objetivos podría ser más razonable la edad al primer hijo.

media_o_na <- function(x) if (all(is.na(x))) NA_real_ else mean(x)

edades_por_familia <- kid_fam[, .(
  n_children          = .N,
  mother_age_at_birth = media_o_na(mother_age_at_birth),
  father_age_at_birth = media_o_na(father_age_at_birth)), by = family_id]

families_edades <- merge(families_raw, edades_por_familia, by = "family_id", all.x = TRUE)
c(familias = nrow(families_raw), filas_resultado = nrow(families_edades))
##        familias filas_resultado 
##           34708           34708

La tabla resultante vuelve a tener como unidad de observación la familia (una fila por familia).

3 Estadística descriptiva y visualización

3.1 Ejercicio 1. Carga y revisión de iris

El archivo usa los nombres SepalLengthCm, SepalWidthCm, PetalLengthCm, PetalWidthCm y Species, equivalentes a Sepal.Length, Sepal.Width, Petal.Length, Petal.Width y Species. Se conservan tal cual. El objeto se llama iris_df para no enmascarar el iris que incluye R.

iris_df <- as.data.frame(iris_data)
str(iris_df)
## 'data.frame':    150 obs. of  6 variables:
##  $ Id           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ SepalLengthCm: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ SepalWidthCm : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ PetalLengthCm: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ PetalWidthCm : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species      : chr  "setosa" "setosa" "setosa" "setosa" ...

Guardar, limpiar y volver a cargar permite comprobar que el objeto se recupera intacto:

rda_file <- tempfile(fileext = ".rda")
save(iris_df, file = rda_file)
rm(iris_df)
load(rda_file)

dim(iris_df)
## [1] 150   6
iris_dt <- as.data.table(iris_df)
class(iris_dt)
## [1] "data.table" "data.frame"

El conjunto tiene 150 observaciones (una flor por fila) y 6 columnas: un identificador, cuatro medidas en cm y la especie.

3.2 Ejercicio 2. Subconjuntos, filtros y variables derivadas

# Solo tres variables
iris_seleccion <- iris_df[, c("SepalLengthCm", "SepalWidthCm", "PetalLengthCm")]
dim(iris_seleccion)
## [1] 150   3
# Todas las variables, solo flores con Sepal.Length >= 5,0
iris_filtrado <- iris_df[iris_df$SepalLengthCm >= 5.0, ]
dim(iris_filtrado)
## [1] 128   6

El filtro conserva 128 de las 150 flores. A continuación se crean las dos variables derivadas con data.table y con dplyr, y se comprueba que dan lo mismo.

# data.table (copy() evita modificar iris_dt por referencia)
iris_dt2 <- copy(iris_dt)
iris_dt2[, `:=`(razon      = SepalLengthCm / SepalWidthCm,
                diferencia = SepalLengthCm - PetalLengthCm)]

# dplyr
iris_dplyr <- iris_df %>%
  mutate(razon      = SepalLengthCm / SepalWidthCm,
         diferencia = SepalLengthCm - PetalLengthCm)

all.equal(as.data.frame(iris_dt2), iris_dplyr)
## [1] TRUE

Estadísticos descriptivos para el conjunto completo y por especie:

resumen <- function(x) {
  c(N = sum(!is.na(x)), media = mean(x, na.rm = TRUE),
    mediana = median(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE))
}
vars_iris <- c("SepalLengthCm", "SepalWidthCm", "PetalLengthCm", "PetalWidthCm",
               "razon", "diferencia")

desc_total <- iris_dt2[, rbindlist(lapply(.SD, function(x) as.list(resumen(x))),
                                   idcol = "variable"), .SDcols = vars_iris]
kable(desc_total, digits = 2, caption = "Descriptivos de iris: conjunto completo")
Descriptivos de iris: conjunto completo
variable N media mediana sd
SepalLengthCm 150 5.84 5.80 0.83
SepalWidthCm 150 3.06 3.00 0.44
PetalLengthCm 150 3.76 4.35 1.77
PetalWidthCm 150 1.20 1.30 0.76
razon 150 1.95 2.03 0.40
diferencia 150 2.09 1.70 1.12 
desc_species <- iris_dt2[, rbindlist(lapply(.SD, function(x) as.list(resumen(x))),
                                     idcol = "variable"),
                         by = Species, .SDcols = vars_iris]
kable(desc_species, digits = 2, caption = "Descriptivos de iris por especie")
Descriptivos de iris por especie
Species variable N media mediana sd
setosa SepalLengthCm 50 5.01 5.00 0.35
setosa SepalWidthCm 50 3.43 3.40 0.38
setosa PetalLengthCm 50 1.46 1.50 0.17
setosa PetalWidthCm 50 0.25 0.20 0.11
setosa razon 50 1.47 1.46 0.12
setosa diferencia 50 3.54 3.50 0.35
versicolor SepalLengthCm 50 5.94 5.90 0.52
versicolor SepalWidthCm 50 2.77 2.80 0.31
versicolor PetalLengthCm 50 4.26 4.35 0.47
versicolor PetalWidthCm 50 1.33 1.30 0.20
versicolor razon 50 2.16 2.16 0.23
versicolor diferencia 50 1.68 1.70 0.35
virginica SepalLengthCm 50 6.59 6.50 0.64
virginica SepalWidthCm 50 2.97 3.00 0.32
virginica PetalLengthCm 50 5.55 5.55 0.55
virginica PetalWidthCm 50 2.03 2.00 0.27
virginica razon 50 2.23 2.17 0.25
virginica diferencia 50 1.04 1.00 0.32

Los promedios globales ocultan diferencias claras entre especies: por ejemplo, el largo del pétalo es muy distinto entre grupos, y por eso la desviación típica del conjunto completo es mayor que la de cada especie por separado. Mirar solo el total llevaría a describir una “flor media” que no se parece a ninguna especie concreta.

3.3 Ejercicio 5. business_staff

3.3.1 Descripción de las variables

str(business_staff)
## Classes 'data.table' and 'data.frame':   1000 obs. of  7 variables:
##  $ id         : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ wave       : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ gender     : Factor w/ 2 levels "female","male": 2 1 2 1 1 2 2 2 1 2 ...
##  $ education  : Factor w/ 3 levels "low","middle",..: 2 2 2 1 1 2 2 3 2 1 ...
##  $ income     : int  3249 3669 4088 1260 1809 791 2114 5456 2627 1153 ...
##  $ happiness  : int  0 7 10 5 3 5 5 3 9 6 ...
##  $ trip_abroad: int  0 1 1 0 0 0 0 1 0 0 ...
##  - attr(*, ".internal.selfref")=<pointer: 0x000001e79624dd90>
summary(business_staff[, .(education, income, happiness, gender, trip_abroad)])
##   education       income       happiness         gender     trip_abroad   
##  low   :221   Min.   :   0   Min.   : 0.000   female:488   Min.   :0.000  
##  middle:582   1st Qu.:2019   1st Qu.: 4.000   male  :512   1st Qu.:0.000  
##  high  :197   Median :3142   Median : 5.000                Median :0.000  
##               Mean   :3230   Mean   : 5.482                Mean   :0.251  
##               3rd Qu.:4240   3rd Qu.: 7.000                3rd Qu.:1.000  
##               Max.   :8731   Max.   :10.000                Max.   :1.000
  • education: nivel educativo (categórica ordinal: low, middle, high).
  • income: ingresos mensuales en euros (numérica, métrica).
  • happiness: felicidad autodeclarada en escala de 0 a 10 (entera, ordinal en origen, usada aquí como cuasi-métrica).
  • gender: female o male (categórica).
  • trip_abroad: indicador 0/1 de si la persona ha realizado un viaje al extranjero.
desc_edu <- business_staff[, .(
  N = .N,
  income_media = mean(income), income_sd = sd(income),
  happiness_media = mean(happiness), happiness_sd = sd(happiness),
  pct_viaje = 100 * mean(trip_abroad)), by = education][order(education)]
kable(desc_edu, digits = 2, caption = "Descriptivos por nivel educativo")
Descriptivos por nivel educativo
education N income_media income_sd happiness_media happiness_sd pct_viaje
low 221 1620.55 497.80 4.52 2.11 23.53
middle 582 3145.97 1006.94 5.61 2.24 24.91
high 197 5285.62 986.58 6.18 2.21 27.41 
kable(as.data.frame(table(gender = business_staff$gender)), caption = "Distribución por género")
Distribución por género
gender Freq
female 488
male 512

3.3.2 Pregunta

¿Cómo se relaciona la felicidad con los ingresos y cambia esa relación según el nivel educativo?

ggplot(business_staff, aes(x = income, y = happiness)) +
  geom_jitter(alpha = 0.35, height = 0.15) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Felicidad e ingresos", x = "Ingresos mensuales (euros)",
       y = "Felicidad (0-10)") +
  theme_minimal()

ggplot(business_staff, aes(x = education, y = happiness)) +
  geom_boxplot() +
  labs(title = "Felicidad según nivel educativo", x = "Nivel educativo",
       y = "Felicidad (0-10)") +
  theme_minimal()

Los ingresos medios aumentan claramente con la educación (1621 < 3146 < 5286 euros en low, middle y high), y la felicidad media también (4.52 < 5.61 < 6.18). El primer gráfico muestra una relación positiva con mucha dispersión. El segundo muestra que, sin controlar por ingresos, la felicidad aumenta con el nivel educativo. Más adelante se verá que esta lectura cambia cuando se comparan personas con el mismo ingreso.

3.3.3 Modelos

model1 <- lm(happiness ~ income, data = business_staff)
model2 <- lm(happiness ~ income + education, data = business_staff)
modelA <- lm(happiness ~ education, data = business_staff)
modelB <- lm(happiness ~ education + gender, data = business_staff)
modelC <- lm(happiness ~ education + gender + income, data = business_staff)
modelD <- lm(happiness ~ income * education, data = business_staff)

modelsummary(list("Model 1" = model1, "Model 2" = model2, "Model A" = modelA,
                  "Model B" = modelB, "Model C" = modelC, "Model D" = modelD),
             estimate  = "{estimate}{stars}",
             statistic = "({std.error})",
             stars     = c("*" = 0.05, "**" = 0.01, "***" = 0.001),
             gof_map   = c("nobs", "r.squared", "adj.r.squared", "AIC", "BIC", "rmse"),
             title     = "Comparación de los modelos de regresión lineal")
Comparación de los modelos de regresión lineal
Model 1 Model 2 Model A Model B Model C Model D
(Intercept) 3.501*** 3.052*** 4.520*** 4.612*** 3.112*** 1.639***
(0.156) (0.179) (0.148) (0.163) (0.192) (0.467)
income 0.001*** 0.001*** 0.001*** 0.002***
(0.000) (0.000) (0.000) (0.000)
educationmiddle -0.293 1.090*** 1.099*** -0.282 1.260*
(0.194) (0.174) (0.174) (0.195) (0.543)
educationhigh -1.658*** 1.662*** 1.660*** -1.649*** 0.347
(0.328) (0.216) (0.216) (0.328) (0.920)
gendermale -0.189 -0.114
(0.140) (0.130)
income × educationmiddle -0.001**
(0.000)
income × educationhigh -0.001**
(0.000)
Num.Obs. 1000 1000 1000 1000 1000 1000
R2 0.164 0.193 0.060 0.062 0.194 0.202
R2 Adj. 0.163 0.191 0.058 0.059 0.190 0.198
RMSE 2.08 2.04 2.20 2.20 2.04 2.03
Modelo Variables independientes Característica principal
Modelo 1 income Solo ingresos
Modelo 2 income + education Ingresos y educación
Modelo A education Solo educación
Modelo B education + gender Educación y género
Modelo C education + gender + income Educación, género e ingresos
Modelo D income * education Ingresos, educación e interacción

3.3.3.1 Modelo 1: felicidad ~ ingresos

Cada euro adicional de ingreso mensual se asocia con 6.1^{-4} puntos más de felicidad, es decir, unos 0.61 puntos por cada 1.000 euros. El modelo explica el 16.4 % de la variabilidad de la felicidad, por lo que los ingresos importan pero dejan sin explicar la mayor parte.

Hay varios problemas al usar esta escala en una regresión lineal. La felicidad está acotada entre 0 y 10, pero el modelo puede predecir fuera de ese rango:

predict(model1, newdata = data.frame(income = 1000000))
##        1 
## 616.6562

Con un millón de euros al mes, el modelo predice 617 puntos, un valor imposible. Además, la escala es subjetiva y discreta, y no tiene por qué haber la misma distancia entre un 4 y un 5 que entre un 8 y un 9. Por último, una asociación no demuestra que los ingresos causen mayor felicidad.

3.3.3.2 Modelos A, B y C: incorporación progresiva de variables independientes

Con low como referencia, el Modelo A estima que los niveles middle y high difieren de low en 1.09 y 1.66 puntos respectivamente.

El género (Modelo B) apenas cambia estos coeficientes y su efecto es pequeño e incierto (coeficiente de gendermale = -0.19).

Al añadir los ingresos (Modelo C), el R² pasa de 0.062 a 0.194 y el R² ajustado de 0.059 a 0.19; el AIC y el BIC también mejoran, por lo que los ingresos aportan mucha más información que el género.

Un matiz importante: la asociación de la educación con la felicidad cambia de signo al controlar por ingresos. Sin controlarlos (Modelo A), mayor educación va con mayor felicidad (1.66 puntos de high respecto a low). Con los ingresos en el modelo (Modelo C), esa diferencia pasa a -1.65: a igual ingreso, las personas con más educación declaran menos felicidad. Esto es coherente con que la educación alta va asociada a ingresos mucho más altos (ver tabla descriptiva): buena parte de la ventaja bruta en felicidad de los más educados se debe a sus ingresos. Los coeficientes de educación en el Modelo C son diferencias condicionadas al ingreso y al género, y no deben leerse como el efecto total de la educación.

¿Debe el modelo final incluir las tres variables? Atendiendo a los criterios, el género no mejora el ajuste ni es relevante, mientras que educación e ingresos sí; el modelo preferible entre A, B y C en términos parsimoniosos sería el que contiene education e income (Modelo 2). El Modelo C solo se justificaría si el género interesara por razones teóricas.

3.3.3.3 Modelo D: interacción ingresos × educación

Qué pasaría si considerasemos que los ingresos estan influidos por la educación, y afectan a la felicidad.

b <- coef(modelD)
pendientes <- c(low    = unname(b["income"]),
                middle = unname(b["income"] + b["income:educationmiddle"]),
                high   = unname(b["income"] + b["income:educationhigh"]))
kable(data.frame(educacion = names(pendientes),
                 pendiente_por_1000_euros = round(pendientes * 1000, 2)),
      caption = "Pendiente estimada de la felicidad por cada 1.000 euros (Modelo D)")
Pendiente estimada de la felicidad por cada 1.000 euros (Modelo D)
educacion pendiente_por_1000_euros
low low 1.78
middle middle 0.86
high high 0.79 
anova(model2, modelD)
## Analysis of Variance Table
## 
## Model 1: happiness ~ income + education
## Model 2: happiness ~ income * education
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)   
## 1    996 4162.1                              
## 2    994 4117.1  2    45.008 5.4332 0.0045 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El efecto de los ingresos es distinto según el nivel educativo: la pendiente es mayor en low (1.78 puntos por 1.000 euros) que en middle (0.86) y high (0.79). El contraste ANOVA entre el Modelo 2 y el Modelo D da F = 5.43 con p = 0.0045, de modo que la interacción mejora significativamente el ajuste. El R² sube de 0.193 a 0.202, una mejora real pero modesta.

rejilla <- expand.grid(
  income = seq(min(business_staff$income), max(business_staff$income), length.out = 100),
  education = levels(business_staff$education))
pred_graf <- rbind(
  transform(rejilla, modelo = "Modelo 2 (sin interacción)", happiness = predict(model2, rejilla)),
  transform(rejilla, modelo = "Modelo D (con interacción)", happiness = predict(modelD, rejilla)))

ggplot(pred_graf, aes(x = income, y = happiness, colour = education)) +
  geom_line(linewidth = 1.1) +
  facet_wrap(~ modelo) +
  labs(title = "Felicidad estimada según ingresos y educación",
       x = "Ingresos mensuales (euros)", y = "Felicidad estimada",
       colour = "Nivel educativo") +
  theme_minimal()

Sin interacción las rectas son paralelas; con interacción la recta de low es más empinada. Modelo elegido: el Modelo D, porque el contraste formal respalda la interacción y mejora el ajuste. Aun así, la mejora es moderada y el modelo es más complejo (y el BIC, más exigente con la complejidad, penaliza su mayor número de parámetros), así que el Modelo 2 sigue siendo una opción más sencilla razonable si interesa la parsimonia.

4 Inferencia básica

4.1 Ejercicio 6. Correlación de Pearson entre income y happiness

cor_test <- cor.test(business_staff$income, business_staff$happiness, method = "pearson")
cor_test
## 
##  Pearson's product-moment correlation
## 
## data:  business_staff$income and business_staff$happiness
## t = 13.968, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3511973 0.4549734
## sample estimates:
##       cor 
## 0.4043861

La correlación es positiva (r = 0.404) y de magnitud moderada: las personas con más ingresos tienden a declarar mayor felicidad. El intervalo de confianza al 95 % va de 0.35 a 0.45, lo que indica que la asociación es claramente positiva y no pequeña. Su cuadrado, r² = 0.164, indica que ingresos y felicidad comparten aproximadamente el 16 % de su varianza, y el resto queda sin explicar.

El p-valor (<2e-16) indica que sería extremadamente raro observar una correlación así si en la población no hubiera relación lineal. Pero eso solo habla de la incertidumbre muestral, no de la importancia práctica ni de la causalidad: lo informativo es el tamaño de r y su intervalo, no que p cruce un umbral. Además, Pearson mide solo la relación lineal y es sensible a valores extremos.

4.2 Ejercicio 7. Prueba t y chi-cuadrado

4.2.1 Ingresos por género

t_test <- t.test(income ~ gender, data = business_staff)
t_test
## 
##  Welch Two Sample t-test
## 
## data:  income by gender
## t = 1.1798, df = 988.29, p-value = 0.2384
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -74.26356 298.19958
## sample estimates:
## mean in group female   mean in group male 
##             3287.693             3175.725

R aplica por defecto la prueba t de Welch, que no asume varianzas iguales. Las mujeres ganan de media 3288 euros y los hombres 3176, una diferencia observada de 112 euros a favor de las mujeres. Sin embargo, el intervalo de confianza al 95 % de la diferencia (-74 a 298 euros) incluye el cero, y el p-valor es 0.238. La diferencia es compatible tanto con una brecha en un sentido, en el otro, o con ninguna: la muestra no permite concluir una diferencia clara, pero tampoco demuestra que no exista.

4.2.2 Educación y género

tabla_edu_gen <- table(education = business_staff$education, gender = business_staff$gender)
kable(tabla_edu_gen, caption = "Frecuencias observadas")
Frecuencias observadas
female male
low 114 107
middle 270 312
high 104 93 
chi <- chisq.test(tabla_edu_gen)
chi
## 
##  Pearson's Chi-squared test
## 
## data:  tabla_edu_gen
## X-squared = 3.2928, df = 2, p-value = 0.1927
kable(round(chi$expected, 1), caption = "Frecuencias esperadas bajo independencia")
Frecuencias esperadas bajo independencia
female male
low 107.8 113.2
middle 284.0 298.0
high 96.1 100.9 
kable(round(chi$stdres, 2), caption = "Residuos estandarizados")
Residuos estandarizados
female male
low 0.94 -0.94
middle -1.80 1.80
high 1.25 -1.25 
v_cramer <- sqrt(unname(chi$statistic) / (sum(tabla_edu_gen) * (min(dim(tabla_edu_gen)) - 1)))

Las frecuencias observadas se parecen a las esperadas si no hubiera asociación; los residuos estandarizados son pequeños (ninguno supera claramente ±2). El estadístico es X² = 3.29 con 2 grados de libertad y p = 0.193, y la V de Cramér es de solo 0.057. La asociación observada entre educación y género es, por tanto, muy débil y compatible con el azar muestral. Con estos datos no hay evidencia de que la distribución educativa difiera entre hombres y mujeres, aunque eso no excluye diferencias pequeñas.

5 Reflexión final

5.1 Ejercicio 8. Estructura y calidad de los datos

Estructura. Lo más importante fue decidir qué es una observación en cada tabla (Wickham, 2014). children y kid_fam tienen como unidad al hijo; families a la familia. Al intentar llevar las edades de los hijos a la tabla de familias hubo que agregar antes (media por familia), una decisión que afecta al resultado: la edad media por familia no es lo mismo que la edad media por hijo, porque las familias numerosas pesan menos. También fue clave comprobar que family_id es único en families y que la fusión conserva los 80.486 hijos. Mantener los datos originales intactos, con diccionario y nombres consistentes, sigue las recomendaciones de Broman y Woo (2018).

Error de medida (Innes et al., 2021). Un ejemplo concreto está en las fechas históricas:

fechas_nac <- c(families_raw$husband_birth_date, families_raw$wife_birth_date,
                children_raw$child_birth_date)
fechas_nac <- fechas_nac[!is.na(fechas_nac)]
prop_enero <- mean(month(fechas_nac) == 1 & mday(fechas_nac) == 1)
c(observada = prop_enero, esperada_si_aleatoria = 1 / 365.25)
##             observada esperada_si_aleatoria 
##           0.141253309           0.002737851

El 14.1 % de las fechas de nacimiento caen el 1 de enero, frente al ~0,27 % esperable si los días se repartieran al azar. Esto sugiere que, cuando solo se conocía el año, se imputó el 1 de enero. Es un error de medida en las fechas de nacimiento que añade ruido a mother_age_at_birth y father_age_at_birth: cada edad puede desviarse hasta casi un año, y como no está distribuido por igual entre familias (los registros más antiguos probablemente tienen más imputaciones), puede afectar a comparaciones entre épocas. Las medias globales son poco sensibles, pero la dispersión y las comparaciones finas sí lo son.

En business_staff, la variable happiness es autodeclarada y se mide con una escala discreta de 0 a 10; las personas pueden usar la escala de forma distinta. Un error de medida no sistemático en una variable resultado aumenta el ruido y reduce la capacidad explicativa (parte del bajo R²), mientras que un error en income (la variable explicativa) tendería a atenuar su coeficiente hacia cero. Por eso, las magnitudes estimadas (por ejemplo, el efecto por cada 1.000 euros) deben leerse como aproximaciones y no como efectos exactos.

6 Referencias

7 Apéndice: información de la sesión

sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Europe/Madrid
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.52          modelsummary_2.6.0  lubridate_1.9.5    
## [4] ggplot2_4.0.3       dplyr_1.2.1         data.table_1.18.6.1
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     lattice_0.22-9     digest_0.6.39     
##  [5] magrittr_2.0.5     evaluate_1.0.5     grid_4.6.1         timechange_0.4.0  
##  [9] RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0     Matrix_1.7-5      
## [13] backports_1.5.1    mgcv_1.9-4         scales_1.4.0       jquerylib_0.1.4   
## [17] cli_3.6.6          rlang_1.3.0        performance_0.19.0 splines_4.6.1     
## [21] tinytable_0.19.0   withr_3.0.3        cachem_1.1.0       yaml_2.3.12       
## [25] otel_0.2.0         datawizard_1.4.1   tools_4.6.1        checkmate_2.3.4   
## [29] bayestestR_0.20.0  vctrs_0.7.3        R6_2.6.1           zoo_1.9-1         
## [33] lifecycle_1.0.5    insight_1.7.0      pkgconfig_2.0.3    pillar_1.11.1     
## [37] bslib_0.12.0       gtable_0.3.6       glue_1.8.1         xfun_0.61         
## [41] tibble_3.3.1       lmtest_0.9-40      tidyselect_1.2.1   rstudioapi_0.19.0 
## [45] parameters_0.30.0  farver_2.1.2       htmltools_0.5.9    nlme_3.1-169      
## [49] tables_0.9.35      rmarkdown_2.32     labeling_0.4.3     compiler_4.6.1    
## [53] S7_0.2.2