Actividad 1

Instalar paquetes y llamar librerías

# ACTIVIDAD EXTRA: app interactiva de datos panel publicada en
# https://warm-florentine-cc0352.netlify.app/

#install.packages("readxl")
library(readxl)
#install.packages("plm")
library(plm)
#install.packages("dplyr")
library(dplyr)
#install.packages("tidyr")
library(tidyr)
#install.packages("jsonlite")
library(jsonlite)

Parte 1. Patentes

Patentes solicitadas por 226 empresas entre 2012 y 2021.

INSTRUCCIONES: Genera el mejor modelo de la Base de Datos “Patentes” y genera predicciones. Incluye conclusiones.

Esta base es la del artículo de Hausman, Hall y Griliches (1984), “Econometric Models for Count Data with an Application to the Patents-R&D Relationship”, Econometrica 52(4), 909-938, así que el modelo se arma siguiendo de cerca lo que ellos hicieron, sin meterse a los modelos de conteo que ahí proponen porque eso ya se sale del módulo. De su trabajo se toman tres decisiones:

1. Un solo control de tamaño y no tres. Ellos reportan que al agregar variables propias de la empresa se les va casi toda la correlación positiva entre I+D y patentes. Aquí pasa lo mismo porque rnd, employ y sales miden todas lo grande que es la empresa: están correlacionadas entre 0.79 y 0.86. Se deja sales como único control de escala y se sacan employ y rnd en nivel.
2. Stock de I+D en lugar del gasto del año. Las patentes salen del conocimiento acumulado y no de lo que se gastó en enero, y esa es la razón por la que ellos trabajan con rezagos. La base ya trae rndstck, que es justo eso.
3. Una tendencia de tiempo. El hallazgo nuevo del artículo es una tendencia negativa en la relación patentes-I+D, o sea que cada año se sacan menos patentes por el mismo dinero invertido. Si eso es cierto, la tendencia tiene que aparecer en el modelo, así que se agrega.

Además todo entra en logaritmos, que es la forma en la que ellos escriben el modelo, y así el coeficiente se lee como elasticidad. También se quita 2021 porque viene incompleto.

pat <- read_excel("PATENT 3.xlsx")
pat <- as.data.frame(pat)

# 2021 viene incompleto: la media de patentes cae de 18.4 a 4.9 y el 55% queda en cero
pat <- pat[pat$year <= 2020, ]
pat <- na.omit(pat[,c("cusip","year","patents","rnd","rndstck","employ","sales","merger")])

# Las tres variables de tamaño miden lo mismo
cor(pat[,c("rnd","employ","sales","rndstck")])

##               rnd    employ     sales   rndstck
## rnd     1.0000000 0.8648309 0.8101436 0.9917219
## employ  0.8648309 1.0000000 0.8340235 0.8582071
## sales   0.8101436 0.8340235 1.0000000 0.8014030
## rndstck 0.9917219 0.8582071 0.8014030 1.0000000

pat$lpat <- log(pat$patents + 1)
pat$lstk <- log(pat$rndstck + 1)
pat$lsal <- log(pat$sales + 1)
pat$t    <- pat$year - 2012

pat <- pdata.frame(pat, index = c("cusip","year"))

# Opción 1 - Modelo de Regresión Agrupada (Pooled)
pooled <- plm(lpat ~ lstk + lsal + merger + t, data=pat, model="pooling")
summary(pooled)

## Pooling Model
## 
## Call:
## plm(formula = lpat ~ lstk + lsal + merger + t, data = pat, model = "pooling")
## 
## Unbalanced Panel: n = 215, T = 2-9, N = 1877
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -2.7288 -0.5513  0.0267  0.6126  2.4076 
## 
## Coefficients:
##              Estimate Std. Error  t-value  Pr(>|t|)    
## (Intercept) -0.127249   0.071111  -1.7894 0.0737063 .  
## lstk         0.675831   0.022171  30.4822 < 2.2e-16 ***
## lsal         0.076663   0.021706   3.5319 0.0004226 ***
## merger       0.491538   0.146644   3.3519 0.0008185 ***
## t           -0.111705   0.007644 -14.6133 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4569.9
## Residual Sum of Squares: 1293.8
## R-Squared:      0.71689
## Adj. R-Squared: 0.71629
## F-statistic: 1185.09 on 4 and 1872 DF, p-value: < 2.22e-16

# Opción 2 - Modelo de Efectos Fijos (Within)
within <- plm(lpat ~ lstk + lsal + merger + t, data=pat, model="within")
summary(within)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = lpat ~ lstk + lsal + merger + t, data = pat, model = "within")
## 
## Unbalanced Panel: n = 215, T = 2-9, N = 1877
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.98282 -0.27653 -0.00397  0.27385  1.44603 
## 
## Coefficients:
##         Estimate Std. Error t-value  Pr(>|t|)    
## lstk    0.184956   0.077950  2.3727 0.0177701 *  
## lsal    0.185259   0.055586  3.3329 0.0008785 ***
## merger  0.040587   0.113024  0.3591 0.7195690    
## t      -0.078950   0.009783 -8.0701 1.336e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    404.08
## Residual Sum of Squares: 380.97
## R-Squared:      0.057198
## Adj. R-Squared: -0.066765
## F-statistic: 25.1471 on 4 and 1658 DF, p-value: < 2.22e-16

# Prueba F
# INTERPRETACIÓN: Si p<0.05 No usar POOLED. Si p>0.05 usar Pooled.
pFtest(within,pooled) #Ojo: El primer argumento es el modelo within!

## 
##  F test for individual effects
## 
## data:  lpat ~ lstk + lsal + merger + t
## F = 18.563, df1 = 214, df2 = 1658, p-value < 2.2e-16
## alternative hypothesis: significant effects

# Opción 3 - Modelo de Efectos Aleatorios (Random)
random <- plm(lpat ~ lstk + lsal + merger + t, data=pat, model="random")
summary(random)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = lpat ~ lstk + lsal + merger + t, data = pat, model = "random")
## 
## Unbalanced Panel: n = 215, T = 2-9, N = 1877
## 
## Effects:
##                  var std.dev share
## idiosyncratic 0.2298  0.4794 0.333
## individual    0.4607  0.6788 0.667
## theta:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.5532  0.7709  0.7709  0.7681  0.7709  0.7709 
## 
## Residuals:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.98e+00 -3.14e-01  1.71e-02 -3.84e-05  3.16e-01  1.45e+00 
## 
## Coefficients:
##              Estimate Std. Error  z-value  Pr(>|z|)    
## (Intercept) -0.250649   0.132823  -1.8871   0.05915 .  
## lstk         0.541346   0.040730  13.2912 < 2.2e-16 ***
## lsal         0.179148   0.038107   4.7012 2.587e-06 ***
## merger       0.083244   0.111408   0.7472   0.45494    
## t           -0.113652   0.005314 -21.3871 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    624.5
## Residual Sum of Squares: 436.05
## R-Squared:      0.30176
## Adj. R-Squared: 0.30027
## Chisq: 824.39 on 4 DF, p-value: < 2.22e-16

# Prueba de Hausman
# INTERPRETACIÓN: Si p<0.05 usar Efectos Fijos. Si p>0.05 usar Efectos Aleatorios.
phtest(random,within) #Ojo: El primer argumento es el modelo random!

## 
##  Hausman Test
## 
## data:  lpat ~ lstk + lsal + merger + t
## chisq = 33.921, df = 4, p-value = 7.734e-07
## alternative hypothesis: one model is inconsistent

# Por lo tanto, el mejor modelo para este panel es el de EFECTOS FIJOS.

Comparación contra el modelo sin ajustar:

viejo <- plm(patents ~ rnd + employ + sales + merger, data=pat, model="within")

data.frame(
  Modelo      = c("Original (niveles, 3 controles de tamaño)",
                  "Ajustado al artículo (logs, stock de I+D, tendencia)"),
  Coef_IyD    = c(round(coef(viejo)["rnd"], 4), round(coef(within)["lstk"], 4)),
  Signo       = c(ifelse(coef(viejo)["rnd"] > 0, "positivo", "negativo"),
                  ifelse(coef(within)["lstk"] > 0, "positivo", "negativo"))
)

##                                                    Modelo Coef_IyD    Signo
## rnd             Original (niveles, 3 controles de tamaño)  -0.0595 negativo
## lstk Ajustado al artículo (logs, stock de I+D, tendencia)   0.1850 positivo

Pronóstico: ¿cuántas patentes solicitarían las 5 empresas más grandes si suben 10% su gasto en I+D?

pat2 <- as.data.frame(pat)
pat2$cusip <- as.character(pat2$cusip)
pat2$year  <- as.numeric(as.character(pat2$year))

top5 <- names(sort(tapply(pat2$patents, pat2$cusip, mean), decreasing=TRUE))[1:5]
pat_pronostico <- pat2[pat2$cusip %in% top5 & pat2$year==2020, ]

intercepto <- fixef(within)
pendiente  <- coef(within)

# El modelo está en logaritmos, así que hay que regresar la predicción a patentes
predecir <- function(stock) {
  l <- intercepto[pat_pronostico$cusip] +
       pendiente["lstk"]*log(stock + 1) +
       pendiente["lsal"]*pat_pronostico$lsal +
       pendiente["merger"]*pat_pronostico$merger +
       pendiente["t"]*pat_pronostico$t
  exp(l) - 1
}

pat_pronostico$Ajustado   <- round(predecir(pat_pronostico$rndstck), 1)
pat_pronostico$ConIyD10   <- round(predecir(pat_pronostico$rndstck * 1.10), 1)
pat_pronostico$Diferencia <- round(pat_pronostico$ConIyD10 - pat_pronostico$Ajustado, 1)

pat_pronostico[,c("cusip","patents","Ajustado","ConIyD10","Diferencia")]

##       cusip patents Ajustado ConIyD10 Diferencia
## 192  149123     140    218.4    222.3        3.9
## 345  277461     214    225.6    229.6        4.0
## 517  369604     667    708.8    721.4       12.6
## 807  459200     333    398.3    405.4        7.1
## 1214 620076     130    151.2    153.9        2.7

CONCLUSIÓN: El mejor modelo sigue siendo el de Efectos Fijos, porque la prueba F descarta el agrupado y la de Hausman descarta el de efectos aleatorios con un p de 7.7e-07.

Lo que concluye el artículo son dos cosas. La primera es que la relación entre I+D y patentes se ve fuerte cuando se comparan empresas entre sí, pero se desarma cuando se controla por las características propias de cada empresa, o sea que buena parte de esa correlación no es el efecto del gasto sino que las empresas grandes gastan más y patentan más al mismo tiempo. La segunda, que ellos marcan como su hallazgo nuevo, es que hay una tendencia negativa en la relación, o sea que cada año se sacan menos patentes por el mismo dinero invertido, y a eso le llaman caída en la productividad de la I+D.

De ahí salieron las decisiones del modelo. Como el efecto se desarma al controlar por la empresa, no tiene caso meter tres variables que miden lo mismo, así que se quedó una sola de tamaño y el stock acumulado de I+D en lugar del gasto de un año suelto. Y como ellos encuentran una tendencia negativa, se metió una tendencia para ver si aquí también aparecía. Todo en logaritmos, que es como lo escriben ellos.

A lo que llegamos nosotros es que las dos conclusiones del artículo se repiten con datos completamente distintos. Ellos trabajaron con 128 empresas entre 1968 y 1974, y esta base es de 226 empresas entre 2012 y 2020, casi cuarenta años después. Aun así la elasticidad del I+D pasa de 0.676 cuando se comparan empresas entre sí a 0.185 cuando se compara cada empresa consigo misma, o sea que se cae un 73%, que es justo lo que ellos describen. Y la tendencia sale en -0.079 y muy significativa, o sea que cada año se sacan alrededor de 7.6% menos patentes con el mismo stock acumulado, igual que su caída de productividad.

Lo que sacamos por nuestra cuenta es que el coeficiente negativo del primer intento no era un problema de los datos sino de cómo estaba armado el modelo. Con las tres variables de tamaño juntas y en niveles, el gasto en I+D salía en -0.0595 y el pronóstico decía que invertir más daba menos patentes. Armándolo como el artículo el coeficiente se voltea a +0.185 y sale significativo, y ahora un 10% más de stock de I+D da 1.78% más patentes, que es poco pero apunta hacia donde tiene sentido. Las cinco empresas más grandes también salen con más patentes al subir la inversión y no con menos.

Dos advertencias para no venderlo de más. El R² dentro bajó de 0.11 a 0.057, pero esos dos números no se comparan porque la variable dependiente ya no es la misma, antes eran patentes y ahora es su logaritmo. Y la conclusión de método del artículo, que es que las patentes son un conteo y piden Poisson o binomial negativa, esa no la seguimos, porque el módulo es de datos en panel con plm. El modelo lineal en logaritmos es lo más cerca que se puede quedar de su especificación con las herramientas de la clase.

Parte 2. Cuidado de la Piel

Valor del mercado de Beauty and Personal Care en México, 2011-2025.

INSTRUCCIONES: Genera el mejor modelo de la Base de Datos “Market sizes” y genera predicciones. Si tuvieras que invertir en alguna sub-categoría, ¿en cuál lo harías?

market <- read_excel("Market sizes.xlsx", skip = 5)

subcats <- c("Bath and Shower","Deodorants","Depilatories","Fragrances",
             "Hair Care","Men's Grooming","Skin Care","Sun Care")

panel <- market %>%
  filter(Category %in% subcats) %>%
  select(Category, matches("^[0-9]{4}$")) %>%
  mutate(across(everything(), as.character)) %>%   #Ojo: los "-" del Excel hacen que unos años se lean como texto
  pivot_longer(-Category, names_to="Año", values_to="Valor") %>%
  mutate(Año  = as.numeric(Año),
         Valor = as.numeric(na_if(Valor,"-"))) %>%
  filter(!is.na(Valor)) %>%
  mutate(Tendencia = Año - 2010,     # 1 = 2011, 15 = 2025, 16 = 2026
         logValor  = log(Valor))      # con log el coeficiente es la tasa de crecimiento

panel <- pdata.frame(panel, index = c("Category","Año"))

# Opción 1 - Modelo de Regresión Agrupada (Pooled)
pooled2 <- plm(logValor ~ Tendencia, data=panel, model="pooling")
summary(pooled2)

## Pooling Model
## 
## Call:
## plm(formula = logValor ~ Tendencia, data = panel, model = "pooling")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -2.386  -0.416   0.395   0.937   1.282 
## 
## Coefficients:
##             Estimate Std. Error t-value Pr(>|t|)    
## (Intercept) 8.910147   0.238229 37.4017  < 2e-16 ***
## Tendencia   0.067392   0.026202  2.5721  0.01135 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    191.64
## Residual Sum of Squares: 181.46
## R-Squared:      0.053087
## Adj. R-Squared: 0.045063
## F-statistic: 6.6155 on 1 and 118 DF, p-value: 0.01135

# Opción 2 - Modelo de Efectos Fijos (Within)
within2 <- plm(logValor ~ Tendencia, data=panel, model="within")
summary(within2)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = logValor ~ Tendencia, data = panel, model = "within")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.25125 -0.04400  0.00802  0.04081  0.24915 
## 
## Coefficients:
##            Estimate Std. Error t-value  Pr(>|t|)    
## Tendencia 0.0673923  0.0018042  37.353 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    10.983
## Residual Sum of Squares: 0.80937
## R-Squared:      0.92631
## Adj. R-Squared: 0.92099
## F-statistic: 1395.23 on 1 and 111 DF, p-value: < 2.22e-16

# Prueba F
# INTERPRETACIÓN: Si p<0.05 No usar POOLED. Si p>0.05 usar Pooled.
pFtest(within2,pooled2) #Ojo: El primer argumento es el modelo within!

## 
##  F test for individual effects
## 
## data:  logValor ~ Tendencia
## F = 3539.4, df1 = 7, df2 = 111, p-value < 2.2e-16
## alternative hypothesis: significant effects

# Opción 3 - Modelo de Efectos Aleatorios (Random)
random2 <- plm(logValor ~ Tendencia, data=panel, model="random")
summary(random2)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = logValor ~ Tendencia, data = panel, model = "random")
## 
## Balanced Panel: n = 8, T = 15, N = 120
## 
## Effects:
##                    var  std.dev share
## idiosyncratic 0.007292 0.085391 0.004
## individual    1.720022 1.311496 0.996
## theta: 0.9832
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.23757 -0.03954  0.00409  0.04301  0.22468 
## 
## Coefficients:
##              Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 8.9101472  0.4639740  19.204 < 2.2e-16 ***
## Tendencia   0.0673923  0.0018042  37.353 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    11.034
## Residual Sum of Squares: 0.86041
## R-Squared:      0.92202
## Adj. R-Squared: 0.92136
## Chisq: 1395.23 on 1 DF, p-value: < 2.22e-16

# Prueba de Hausman
# INTERPRETACIÓN: Si p<0.05 usar Efectos Fijos. Si p>0.05 usar Efectos Aleatorios.
phtest(random2,within2) #Ojo: El primer argumento es el modelo random!

## 
##  Hausman Test
## 
## data:  logValor ~ Tendencia
## chisq = 1.6241e-14, df = 1, p-value = 1
## alternative hypothesis: one model is inconsistent

# Por lo tanto, el mejor modelo para este panel es el de EFECTOS ALEATORIOS.

Pronóstico 2026 por subcategoría:

panel2 <- as.data.frame(panel)

pronostico <- panel2 %>%
  group_by(Category) %>%
  summarise(Crecimiento = coef(lm(logValor ~ Tendencia))[2],
            R2          = summary(lm(logValor ~ Tendencia))$r.squared,
            Pron2026    = exp(predict(lm(logValor ~ Tendencia),
                                      data.frame(Tendencia=16)))) %>%
  arrange(desc(Pron2026))

as.data.frame(pronostico)

##          Category Crecimiento        R2  Pron2026
## 1       Skin Care  0.07443372 0.9703365 72529.360
## 2       Hair Care  0.05418877 0.9728327 55317.968
## 3      Fragrances  0.07462101 0.9032801 52078.190
## 4  Men's Grooming  0.06504529 0.9616394 50701.874
## 5      Deodorants  0.05670968 0.9230642 22763.378
## 6 Bath and Shower  0.06533531 0.9865601 22102.258
## 7        Sun Care  0.08835282 0.9349442  4412.509
## 8    Depilatories  0.06045172 0.9689007  2147.434

# La sub-categoría que más crece y la que menos crece
pronostico[which.max(pronostico$Crecimiento), ]

## # A tibble: 1 × 4
##   Category Crecimiento    R2 Pron2026
##   <fct>          <dbl> <dbl>    <dbl>
## 1 Sun Care      0.0884 0.935    4413.

pronostico[which.min(pronostico$Crecimiento), ]

## # A tibble: 1 × 4
##   Category  Crecimiento    R2 Pron2026
##   <fct>           <dbl> <dbl>    <dbl>
## 1 Hair Care      0.0542 0.973   55318.

CONCLUSIÓN: Definitivamente invertir en Skin Care. Con el modelo de Efectos Aleatorios, el mercado de BPC en México crece 6.74% anual en promedio, pero no todas las subcategorías van al mismo ritmo: Sun Care crece más rápido (8.84%) y Hair Care más lento (5.42%). Ninguna de las dos es la respuesta obvia: lo que crece mucho atrae más competencia, y lo que crece poco suele ser porque ya está maduro y no le interesa a nadie más. Sun Care además es diminuta (4,413 millones contra 55,318 de Hair Care), así que aunque crezca rápido, el tamaño no compensa. Skin Care es el punto medio que sí conviene: crece casi tan rápido como Sun Care (7.44%), ya es la categoría más grande de las ocho (72,529 millones proyectados para 2026), y es la que mejor le ajusta al modelo (R² de 0.9703), así que es en la que más confío. Eso sí, esta base solo mide tamaño de mercado, no competencia real, así que para cerrar bien el argumento haría falta cruzarlo con la base de Company Shares.

Parte 3. Banco Mundial

Esperanza de vida en 12 países, 2000-2023.

INSTRUCCIONES: Importa datos del banco mundial y genera el mejor modelo, incluyendo predicciones. Incluye conclusiones y gráficas.

options(timeout = 600)

paises <- "MX;US;CA;BR;CL;CO;AR;PE;ES;KR;JP;IN"

baja <- function(indicador) {
  u <- paste0("https://api.worldbank.org/v2/country/", paises,
              "/indicator/", indicador, "?format=json&per_page=2000&date=2000:2023")
  d <- fromJSON(u)[[2]]
  data.frame(Pais = d$country$value, Año = as.numeric(d$date), Valor = d$value)
}

if (file.exists("wdi.rds")) {
  wdi <- readRDS("wdi.rds")
} else {
  esperanza  <- baja("SP.DYN.LE00.IN")     # esperanza de vida al nacer
  pib        <- baja("NY.GDP.PCAP.CD")     # PIB per cápita
  salud      <- baja("SH.XPD.CHEX.GD.ZS")  # gasto en salud % del PIB
  urbana     <- baja("SP.URB.TOTL.IN.ZS")  # población urbana %
  fertilidad <- baja("SP.DYN.TFRT.IN")     # hijos por mujer

  wdi <- esperanza %>%
    rename(Esperanza = Valor) %>%
    left_join(rename(pib,        PIB        = Valor), by=c("Pais","Año")) %>%
    left_join(rename(salud,      Salud      = Valor), by=c("Pais","Año")) %>%
    left_join(rename(urbana,     Urbana     = Valor), by=c("Pais","Año")) %>%
    left_join(rename(fertilidad, Fertilidad = Valor), by=c("Pais","Año")) %>%
    mutate(PIB = PIB/1000) %>%          # PIB per cápita en miles de USD
    na.omit()

  saveRDS(wdi, "wdi.rds")
}

head(wdi)

##        Pais  Año Esperanza       PIB    Salud   Urbana Fertilidad
## 1 Argentina 2023    77.395 14.261847 10.26952 92.19451      1.500
## 2 Argentina 2022    75.806 13.962189 10.15950 92.11284      1.482
## 3 Argentina 2021    73.948 10.738018 11.00910 92.02928      1.585
## 4 Argentina 2020    75.878  8.535599 10.67246 91.94390      1.601
## 5 Argentina 2019    76.847  9.955975 10.11911 91.85678      1.882
## 6 Argentina 2018    76.770 11.752800 10.17289 91.76799      2.067

# Ingreso contra esperanza de vida. La recta roja es el ajuste lineal y la curva azul
# sigue los datos sin imponerle forma, para ver que la relación no es una línea.
plot(wdi$PIB, wdi$Esperanza,
     xlab="PIB per cápita (miles de USD)", ylab="Esperanza de vida (años)",
     main="Ingreso y esperanza de vida, 12 países 2000-2023",
     pch=19, col="grey40")
abline(lm(Esperanza ~ PIB, data=wdi), col="red", lwd=2)
lines(lowess(wdi$PIB, wdi$Esperanza), col="blue", lwd=2)
legend("bottomright", c("Ajuste lineal","Tendencia real"),
       col=c("red","blue"), lwd=2, bty="n")

wdi <- pdata.frame(wdi, index = c("Pais","Año"))

# Opción 1 - Modelo de Regresión Agrupada (Pooled)
pooled3 <- plm(Esperanza ~ PIB + Salud + Urbana + Fertilidad, data=wdi, model="pooling")
summary(pooled3)

## Pooling Model
## 
## Call:
## plm(formula = Esperanza ~ PIB + Salud + Urbana + Fertilidad, 
##     data = wdi, model = "pooling")
## 
## Balanced Panel: n = 12, T = 24, N = 288
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -6.065  -0.595   0.250   0.971   3.527 
## 
## Coefficients:
##               Estimate Std. Error  t-value  Pr(>|t|)    
## (Intercept) 75.4130817  0.9534218  79.0973 < 2.2e-16 ***
## PIB          0.1278894  0.0109406  11.6894 < 2.2e-16 ***
## Salud       -0.3112481  0.0573162  -5.4304 1.212e-07 ***
## Urbana       0.1194947  0.0086605  13.7977 < 2.2e-16 ***
## Fertilidad  -4.3246340  0.2745273 -15.7530 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    6263.9
## Residual Sum of Squares: 811.36
## R-Squared:      0.87047
## Adj. R-Squared: 0.86864
## F-statistic: 475.457 on 4 and 283 DF, p-value: < 2.22e-16

# Opción 2 - Modelo de Efectos Fijos (Within)
within3 <- plm(Esperanza ~ PIB + Salud + Urbana + Fertilidad, data=wdi, model="within")
summary(within3)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = Esperanza ~ PIB + Salud + Urbana + Fertilidad, 
##     data = wdi, model = "within")
## 
## Balanced Panel: n = 12, T = 24, N = 288
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -5.187  -0.387   0.221   0.673   1.981 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## PIB         0.068722   0.014465  4.7508 3.284e-06 ***
## Salud      -0.045645   0.084567 -0.5397    0.5898    
## Urbana      0.338864   0.042202  8.0296 2.967e-14 ***
## Fertilidad -2.302993   0.337533 -6.8230 5.776e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    795.01
## Residual Sum of Squares: 332.54
## R-Squared:      0.58171
## Adj. R-Squared: 0.55864
## F-statistic: 94.5665 on 4 and 272 DF, p-value: < 2.22e-16

# Prueba F
# INTERPRETACIÓN: Si p<0.05 No usar POOLED. Si p>0.05 usar Pooled.
pFtest(within3,pooled3) #Ojo: El primer argumento es el modelo within!

## 
##  F test for individual effects
## 
## data:  Esperanza ~ PIB + Salud + Urbana + Fertilidad
## F = 35.603, df1 = 11, df2 = 272, p-value < 2.2e-16
## alternative hypothesis: significant effects

# Opción 3 - Modelo de Efectos Aleatorios (Random)
random3 <- plm(Esperanza ~ PIB + Salud + Urbana + Fertilidad, data=wdi, model="random")
summary(random3)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = Esperanza ~ PIB + Salud + Urbana + Fertilidad, 
##     data = wdi, model = "random")
## 
## Balanced Panel: n = 12, T = 24, N = 288
## 
## Effects:
##                 var std.dev share
## idiosyncratic 1.223   1.106 0.355
## individual    2.220   1.490 0.645
## theta: 0.8502
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
##  -5.295  -0.450   0.206   0.687   2.515 
## 
## Coefficients:
##              Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 65.275337   2.346653 27.8164 < 2.2e-16 ***
## PIB          0.075145   0.014225  5.2825 1.274e-07 ***
## Salud       -0.022270   0.081097 -0.2746    0.7836    
## Urbana       0.199977   0.027046  7.3939 1.426e-13 ***
## Fertilidad  -2.926314   0.321027 -9.1155 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    917.69
## Residual Sum of Squares: 371.22
## R-Squared:      0.59548
## Adj. R-Squared: 0.58977
## Chisq: 416.602 on 4 DF, p-value: < 2.22e-16

# Prueba de Hausman
# INTERPRETACIÓN: Si p<0.05 usar Efectos Fijos. Si p>0.05 usar Efectos Aleatorios.
phtest(random3,within3) #Ojo: El primer argumento es el modelo random!

## 
##  Hausman Test
## 
## data:  Esperanza ~ PIB + Salud + Urbana + Fertilidad
## chisq = 34.163, df = 4, p-value = 6.9e-07
## alternative hypothesis: one model is inconsistent

# Por lo tanto, el mejor modelo para este panel es el de EFECTOS FIJOS.

Pronóstico: esperanza de vida de cada país si su PIB per cápita sube 10%.

wdi2 <- as.data.frame(wdi)
wdi2$Pais <- as.character(wdi2$Pais)
wdi2$Año <- as.numeric(as.character(wdi2$Año))

wdi_pronostico <- wdi2[wdi2$Año == max(wdi2$Año), ]

intercepto3 <- fixef(within3)
pendiente3  <- coef(within3)

wdi_pronostico$Ajustado <- intercepto3[wdi_pronostico$Pais] +
                           pendiente3["PIB"]*wdi_pronostico$PIB +
                           pendiente3["Salud"]*wdi_pronostico$Salud +
                           pendiente3["Urbana"]*wdi_pronostico$Urbana +
                           pendiente3["Fertilidad"]*wdi_pronostico$Fertilidad

wdi_pronostico$Pronostico <- intercepto3[wdi_pronostico$Pais] +
                             pendiente3["PIB"]*wdi_pronostico$PIB*1.10 +
                             pendiente3["Salud"]*wdi_pronostico$Salud +
                             pendiente3["Urbana"]*wdi_pronostico$Urbana +
                             pendiente3["Fertilidad"]*wdi_pronostico$Fertilidad

wdi_pronostico[,c("Pais","Esperanza","Ajustado","Pronostico")]

##              Pais Esperanza Ajustado Pronostico
## 24      Argentina  77.39500 77.90201   78.00002
## 48         Brazil  75.84800 75.02477   75.09609
## 72         Canada  81.62634 83.00841   83.38533
## 96          Chile  81.16700 81.12941   81.24680
## 120      Colombia  77.72500 75.98704   76.03523
## 144         India  72.00300 70.14654   70.16327
## 168         Japan  84.04122 84.22867   84.47068
## 192   Korea, Rep.  83.42927 81.82440   82.06956
## 216        Mexico  75.06900 75.51749   75.61254
## 240          Peru  77.74000 76.76254   76.81696
## 264         Spain  83.93415 83.17301   83.40318
## 288 United States  78.38537 80.46590   81.03345

CONCLUSIÓN: El mejor modelo es el de Efectos Fijos (la prueba F descarta el agrupado y Hausman da un p de 6.9e-07). Explica el 58.17% de la variación dentro de cada país, y tres de las cuatro variables salen significativas: 1,000 dólares más de PIB per cápita suman 0.069 años de vida, un punto más de población urbana suma 0.339, y un hijo más por mujer resta 2.303. El gasto en salud no salió significativo, lo cual tiene sentido porque gastar más no es lo mismo que gastar bien. Lo que más me llamó la atención es que el ingreso mueve poco la aguja: subir el PIB per cápita 10% le suma a México menos de una décima de año, y al que más le sube, Estados Unidos, le agrega poco más de medio año. La gráfica explica por qué: la relación no es una recta, sino que se aplana después de los primeros miles de dólares. El modelo le atina bien a países como Chile o Japón, pero falla por más de dos años en Estados Unidos, así que sirve para comparar tendencias entre países, no para predecir con precisión la esperanza de vida de uno solo.