Modelización del Radio Planetario en Exoplanetas: Un Enfoque de Regresión Lineal Múltiple

Este proyecto parte de una idea sencilla, pero muy completa desde el punto de vista estadístico: estudiar qué variables físicas, orbitales y estelares están asociadas con el radio de un exoplaneta. En vez de quedarme solo en enseñar unos resultados finales, el objetivo de este artículo es recorrer todo el proceso como un tutorial: cargar la base de datos, seleccionar variables, explorar relaciones, ajustar un modelo lineal múltiple, diagnosticar sus problemas y probar una mejora mediante transformación logarítmica.

La gracia del proyecto está en que reúne muchos conceptos de la asignatura en un caso real: limpieza de datos, análisis exploratorio, correlaciones simples y parciales, multicolinealidad, regresión lineal múltiple, interpretación de coeficientes, intervalos de confianza, bondad del ajuste, análisis de residuos, distancia de Cook, selección automática mediante AIC y transformación logarítmica.

La variable que se quiere explicar es pl_rade, el radio planetario medido en radios terrestres. Es decir, se toma la Tierra como unidad de referencia y se analiza qué factores ayudan a explicar que un exoplaneta sea más pequeño o más grande.

1. Introducción

1.1 Descripción de la base de datos

La base de datos utilizada procede del NASA Exoplanet Archive, concretamente de la tabla Planetary Systems Composite Parameters. Esta base recoge información sobre exoplanetas confirmados y sus sistemas planetarios, incluyendo variables relativas al planeta, a su órbita, a la estrella anfitriona y a la distancia del sistema respecto a la Tierra.

El objetivo del análisis es estudiar qué características físicas, orbitales y estelares están asociadas con el tamaño de un exoplaneta. Para ello, se selecciona como variable dependiente el radio planetario (pl_rade), medido en radios terrestres. Esta variable permite analizar el tamaño relativo de los exoplanetas respecto a la Tierra.

Desde el punto de vista del problema, el modelo puede entenderse como una herramienta exploratoria para estudiar qué factores están relacionados con el tamaño planetario. Este tipo de análisis puede ser útil en el contexto de la investigación astronómica, por ejemplo, para caracterizar poblaciones de exoplanetas, identificar patrones físicos y apoyar la priorización de candidatos para observaciones posteriores.

library(ppcor)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(broom)
library(tibble)
library(lmtest)

exoplanetas <- read_csv(
  "PSCompPars_2026.04.29_08.19.44.csv",
  comment = "#"
)

glimpse(exoplanetas)
summary(exoplanetas)

Rows: 6273 Columns: 84
$ pl_name         <chr> "11 Com b", "11 UMi b", "14 And b", ...
$ hostname        <chr> "11 Com", "11 UMi", "14 And", ...
$ discoverymethod <chr> "Radial Velocity", "Radial Velocity", ...
$ disc_year       <dbl> 2007, 2009, 2008, 2002, 1996, ...
$ pl_orbper       <dbl> 323.21, 516.22, 186.76, 1766.41, ...
$ pl_orbsmax      <dbl> 1.178, 1.530, 0.775, 2.839, ...
$ pl_rade         <dbl> 12.2, 12.3, 13.1, 12.5, ...
$ pl_bmasse       <dbl> 4914.8985, 4684.8142, 1131.1513, ...
$ pl_eqt          <dbl> NA, NA, NA, NA, NA, ...
$ st_teff         <dbl> 4874, 4213, 4888, 5338, ...
$ st_rad          <dbl> 13.76, 29.79, 11.55, 0.93, ...
$ st_mass         <dbl> 2.09, 2.78, 1.78, 0.97, ...
$ sy_dist         <dbl> 93.1846, 125.3210, 75.4392, ...

Interpretación. La base original es amplia: contiene miles de exoplanetas y 84 columnas. No todas sirven directamente para el modelo, porque hay identificadores, errores de medición, límites, unidades alternativas y variables con muchos valores perdidos. Por eso la primera decisión importante no es ajustar el modelo, sino decidir qué variables tienen sentido estadístico y físico.

1.2 Selección y justificación de variables

La base de datos original contiene 84 variables. Sin embargo, no todas se incorporan al modelo de regresión, ya que algunas tienen un papel meramente identificativo, otras representan errores de medición, otras son variables equivalentes expresadas en distintas unidades y otras tienen un elevado número de valores perdidos o una interpretación menos directa para el objetivo del análisis.

El objetivo del estudio es explicar el radio planetario (pl_rade), por lo que se seleccionan variables con una relación física, orbital, estelar u observacional razonable con el tamaño del planeta.

1.2.1 Variables no incorporadas al modelo

En primer lugar, se excluyen las variables identificativas o descriptivas, como pl_name, hostname, disc_facility, rastr y decstr. Estas variables permiten identificar el planeta, la estrella anfitriona, la instalación de descubrimiento o la posición astronómica en formato textual, pero no aportan una magnitud cuantitativa directamente interpretable en un modelo de regresión lineal.

También se excluyen variables que representan la misma magnitud en otra unidad. Por ejemplo, pl_radj mide el radio planetario en radios de Júpiter, mientras que pl_rade lo mide en radios terrestres. Como pl_rade se ha elegido como variable dependiente, incluir pl_radj como explicativa supondría introducir una variable prácticamente equivalente a la respuesta. Del mismo modo, pl_bmassj mide la masa planetaria en masas de Júpiter, mientras que pl_bmasse la mide en masas terrestres; por ello se conserva únicamente pl_bmasse.

Además, se excluyen las variables asociadas a errores de medición, como pl_orbpererr1, pl_orbpererr2, pl_orbsmaxerr1, pl_orbsmaxerr2, pl_radeerr1, pl_radeerr2, pl_bmasseerr1, pl_bmasseerr2, st_tefferr1, st_tefferr2, st_raderr1, st_raderr2, st_masserr1, st_masserr2, sy_disterr1 y sy_disterr2. Estas columnas informan sobre la incertidumbre de las estimaciones, pero no representan directamente características físicas del planeta o de la estrella.

De forma similar, se excluyen las variables terminadas en lim, como pl_orbperlim, pl_orbsmaxlim, pl_radelim, pl_bmasselim, pl_orbeccenlim, pl_eqtlim, st_tefflim, st_radlim, st_masslim, st_metlim y st_logglim. Estas variables indican si una medición corresponde a un límite superior o inferior, por lo que tienen un carácter auxiliar y no se interpretan como variables explicativas principales.

Por último, algunas variables potencialmente interesantes, como pl_orbeccen, pl_insol, st_met, st_logg, sy_vmag, sy_kmag o sy_gaiamag, no se incorporan al modelo inicial para evitar una pérdida excesiva de observaciones por valores perdidos y para mantener un modelo interpretable. No obstante, podrían considerarse en análisis complementarios o ampliaciones posteriores.

1.2.2 Variables seleccionadas para el análisis

Tras la revisión anterior, se selecciona un conjunto reducido de variables con interpretación clara y relación potencial con el radio planetario. La variable dependiente será pl_rade, que representa el radio del planeta en radios terrestres.

datos <- exoplanetas %>%
  dplyr::select(
    disc_year,
    pl_rade,
    pl_bmasse,
    pl_orbper,
    pl_orbsmax,
    pl_eqt,
    st_teff,
    st_rad,
    st_mass,
    sy_dist
  )

datos <- datos %>%
  drop_na() %>%
  filter(
    pl_rade > 0,
    pl_bmasse > 0,
    pl_orbper > 0,
    pl_orbsmax > 0,
    sy_dist > 0
  )

Interpretación. Este bloque deja el conjunto de datos preparado para el análisis. Primero se seleccionan las variables relevantes y después se eliminan filas con valores perdidos. También se filtran valores no positivos en variables donde el valor debe ser positivo para que la interpretación física tenga sentido y para poder aplicar transformaciones logarítmicas más adelante.

2. Análisis Exploratorio de Datos

2.1 Resumen descriptivo

summary(datos)

disc_year     pl_rade        pl_bmasse       pl_orbper
Min.   :1996  Min.   : 0.3098  Min.   :   0.0374  Min.   :        0
1st Qu.:2014  1st Qu.: 1.6400  1st Qu.:   3.5300  1st Qu.:        4
Median :2016  Median : 2.4700  Median :   7.1422  Median :        8
Mean   :2017  Mean   : 4.4750  Mean   : 145.3384  Mean   :    93906
3rd Qu.:2021  3rd Qu.: 4.1975  3rd Qu.:  23.6900  3rd Qu.:       21
Max.   :2026  Max.   :25.0000  Max.   :8899.1954  Max.   :402000000

pl_orbsmax       pl_eqt         st_teff       st_rad
Min.   :   0.0050  Min.   :  55.9  Min.   : 2375  Min.   :0.0131
1st Qu.:   0.0474  1st Qu.: 564.0  1st Qu.: 5018  1st Qu.:0.7700
Median :   0.0780  Median : 812.0  Median : 5598  Median :0.9400
Mean   :   1.9882  Mean   : 907.3  Mean   : 5382  Mean   :1.0083
3rd Qu.:   0.1428  3rd Qu.:1159.8  3rd Qu.: 5931  3rd Qu.:1.1900
Max.   :7506.0000  Max.   :4050.0  Max.   :10170  Max.   :6.3000

st_mass          sy_dist
Min.   :0.0160  Min.   :   1.301
1st Qu.:0.8000  1st Qu.: 202.585
Median :0.9403  Median : 492.485
Mean   :0.9343  Mean   : 589.504
3rd Qu.:1.0800  3rd Qu.: 864.331
Max.   :2.7800  Max.   :3460.510

Interpretación. Este resumen permite detectar rápidamente la escala y la dispersión de cada variable. La mediana del radio planetario es 2.47 radios terrestres, mientras que la media es 4.475. Esa diferencia ya sugiere asimetría positiva: la mayoría de planetas tienen radios relativamente bajos, pero algunos valores grandes empujan la media hacia arriba.

También se aprecia mucha dispersión en pl_bmasse, pl_orbper, pl_orbsmax y sy_dist. Esto indica que la muestra mezcla sistemas planetarios muy distintos, algo esperable en una base astronómica que combina exoplanetas de múltiples tipos.

ggplot(datos, aes(x = pl_rade)) +
  geom_histogram(bins = 40, fill = "#E67E22", color = "white") +
  labs(
    title = "Distribución del radio planetario",
    subtitle = "Variable dependiente del modelo",
    x = "Radio planetario (radios terrestres)",
    y = "Frecuencia"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

La distribución de pl_rade presenta una clara asimetría positiva, con una concentración de planetas de radio reducido y algunos valores extremos de mayor tamaño. Esto anticipa posibles problemas de normalidad en los residuos del modelo lineal.

El análisis descriptivo muestra una elevada dispersión en varias variables. En particular, pl_bmasse, pl_orbper, pl_orbsmax y sy_dist presentan rangos amplios, lo que indica que la muestra incluye exoplanetas y sistemas planetarios muy heterogéneos. También se observa que algunas variables presentan diferencias importantes entre la media y la mediana, lo que sugiere distribuciones asimétricas y posible presencia de valores extremos.

2.2 Análisis gráfico

plot(datos)

Salida. El comando plot(datos) genera una matriz de gráficos de dispersión por pares entre las variables seleccionadas. En el PDF aparece como una matriz grande donde cada panel compara dos variables distintas.

Interpretación. El gráfico de dispersión por pares permite observar visualmente las relaciones bivariantes entre las variables seleccionadas. En general, no se aprecian relaciones lineales fuertes entre pl_rade y la mayoría de variables explicativas. Sin embargo, sí se observa una relación muy marcada entre pl_orbper y pl_orbsmax, lo que anticipa un posible problema de multicolinealidad entre ambas variables orbitales.

2.3 Correlaciones con la variable dependiente

# Correlaciones con la variable dependiente
cor(datos)["pl_rade", ]

disc_year       pl_rade     pl_bmasse     pl_orbper    pl_orbsmax
-0.0008999    1.0000000     0.4240659     0.0264991     0.0283655

pl_eqt        st_teff       st_rad        st_mass       sy_dist
0.4244222     0.2730254     0.4068889     0.4115688    -0.1405485

matriz_cor_simple <- cor(datos, use = "complete.obs")
round(matriz_cor_simple, 3)

disc_year pl_rade pl_bmasse pl_orbper pl_orbsmax pl_eqt st_teff st_rad st_mass sy_dist
disc_year        1.000  -0.001     0.008     0.013      0.013 -0.056  -0.189 -0.086  -0.151  -0.258
pl_rade         -0.001   1.000     0.424     0.026      0.028  0.424   0.273  0.407   0.412  -0.141
pl_bmasse        0.008   0.424     1.000     0.067      0.072  0.184   0.147  0.201   0.208  -0.080
pl_orbper        0.013   0.026     0.067     1.000      1.000 -0.015  -0.037 -0.022  -0.032  -0.019
pl_orbsmax       0.013   0.028     0.072     1.000      1.000 -0.015  -0.036 -0.022  -0.031  -0.019
pl_eqt          -0.056   0.424     0.184    -0.015     -0.015  1.000   0.459  0.463   0.522   0.037
st_teff         -0.189   0.273     0.147    -0.037     -0.036  0.459   1.000  0.683   0.897   0.438
st_rad          -0.086   0.407     0.201    -0.022     -0.022  0.463   0.683  1.000   0.817   0.241
st_mass         -0.151   0.412     0.208    -0.032     -0.031  0.522   0.897  0.817   1.000   0.322
sy_dist         -0.258  -0.141    -0.080    -0.019     -0.019  0.037   0.438  0.241   0.322   1.000

cor_pl_rade <- cor(datos)["pl_rade", ] %>%
  as.data.frame() %>%
  rownames_to_column("variable") %>%
  rename(correlacion = 2) %>%
  filter(variable != "pl_rade") %>%
  mutate(variable = reorder(variable, correlacion))

ggplot(cor_pl_rade, aes(x = correlacion, y = variable, fill = correlacion > 0)) +
  geom_col(show.legend = FALSE) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_fill_manual(values = c("#555555", "#E67E22")) +
  labs(
    title = "Correlación simple con el radio planetario",
    x = "Correlación con pl_rade",
    y = "Variable"
  ) +
  theme_minimal(base_size = 13)

La matriz de correlaciones simples muestra la relación lineal entre cada variable numérica y la variable dependiente pl_rade. Como era esperable, la correlación de pl_rade consigo misma es igual a 1, por lo que no se interpreta.

Las variables que presentan una mayor asociación lineal con el radio planetario son pl_eqt y pl_bmasse, ambas con una correlación aproximada de 0.424. Esto indica una relación positiva moderada: los planetas con mayor temperatura de equilibrio o mayor masa tienden, en promedio, a presentar radios superiores. También se observan correlaciones positivas moderadas con st_mass y st_rad, con valores cercanos a 0.41.

La variable st_teff presenta una correlación positiva más débil, de aproximadamente 0.273, mientras que sy_dist muestra una correlación negativa débil, cercana a -0.141. Por último, pl_orbper, pl_orbsmax y disc_year presentan correlaciones prácticamente nulas con pl_rade. No obstante, pl_orbper y pl_orbsmax se analizarán con más detalle porque pueden presentar multicolinealidad entre ellas.

2.4 Análisis previo de multicolinealidad

Antes de calcular la matriz de correlaciones parciales, se revisa la matriz de correlaciones simples entre las variables explicativas. Este paso es importante porque las correlaciones parciales requieren invertir una matriz de covarianzas/correlaciones, y si existen variables casi perfectamente correlacionadas entre sí pueden aparecer problemas numéricos.

cor(datos$pl_orbper, datos$pl_orbsmax)

[1] 0.9997907

ggplot(datos, aes(x = pl_orbper, y = pl_orbsmax)) +
  geom_point(alpha = 0.35, color = "#E67E22") +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    title = "Relación entre periodo orbital y semieje mayor",
    subtitle = "Escala logarítmica",
    x = "Periodo orbital (días)",
    y = "Semieje mayor orbital (UA)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

El resultado muestra una correlación extremadamente elevada entre pl_orbper y pl_orbsmax, cercana a 1. Esto indica que ambas variables aportan información muy similar. Desde un punto de vista físico, este resultado es coherente, ya que los planetas más alejados de su estrella suelen tener periodos orbitales más largos.

# Intento inicial de calcular correlaciones parciales con todas las variables
pcor(datos)$estimate

Error in solve.default(cvx): sistema es computacionalmente singular: número de condición recíproco = 1.78325e-16

El cálculo anterior puede generar un problema numérico debido a la elevada multicolinealidad entre pl_orbper y pl_orbsmax. En concreto, la función pcor() necesita invertir internamente una matriz, y esta operación puede fallar cuando dos variables están casi perfectamente correlacionadas.

# Comparación de la relación de cada variable orbital con la variable dependiente
cor(datos$pl_rade, datos$pl_orbper)
cor(datos$pl_rade, datos$pl_orbsmax)

[1] 0.02649906
[1] 0.0283655

Como pl_orbsmax presenta una correlación ligeramente superior con pl_rade, se conserva esta variable y se elimina pl_orbper para evitar problemas de multicolinealidad.

datos <- datos %>%
  dplyr::select(
    pl_rade,
    pl_bmasse,
    pl_orbsmax,
    pl_eqt,
    st_teff,
    st_rad,
    st_mass,
    sy_dist
  )

2.5 Correlaciones parciales con la variable dependiente

# Correlaciones parciales con la variable dependiente
pcor(datos)$estimate["pl_rade", ]

pl_rade    pl_bmasse  pl_orbsmax  pl_eqt     st_teff    st_rad     st_mass    sy_dist
1.0000000  0.3408100  0.0177740   0.2269549 -0.1317693 0.0796841 0.1825776 -0.1958625

matriz_cor_parcial <- pcor(datos)$estimate
round(matriz_cor_parcial, 3)

pl_rade pl_bmasse pl_orbsmax pl_eqt st_teff st_rad st_mass sy_dist
pl_rade        1.000     0.341      0.018  0.227  -0.132  0.080   0.183  -0.196
pl_bmasse      0.341     1.000      0.019  0.036   0.020  0.044   0.046  -0.026
pl_orbsmax     0.018     0.019      1.000 -0.006  -0.003 -0.002  -0.006  -0.004
pl_eqt         0.227     0.036     -0.006  1.000   0.027  0.024   0.112  -0.105
st_teff       -0.132     0.020     -0.003  0.027   1.000  0.174   0.791   0.287
st_rad         0.080     0.044     -0.002  0.024   0.174  1.000   0.526  -0.094
st_mass        0.183     0.046     -0.006  0.112   0.791  0.526   1.000  -0.012
sy_dist       -0.196    -0.026     -0.004 -0.105   0.287 -0.094  -0.012   1.000

Las correlaciones parciales muestran que, al controlar por el resto de variables, pl_bmasse sigue siendo la variable más relacionada con pl_rade, con una correlación parcial positiva de 0.341. También mantienen cierta asociación pl_eqt, sy_dist y st_mass, aunque con valores más moderados. En cambio, pl_orbsmax y st_rad presentan correlaciones parciales muy bajas, lo que sugiere que su aportación específica al explicar el radio planetario es limitada.

Además, st_teff cambia de signo respecto a la correlación simple, pasando de una relación positiva a una relación parcial negativa. Esto sugiere que la asociación simple entre la temperatura de la estrella y el radio planetario podía estar condicionada por su relación con otras variables del modelo.

2.6 Conclusión del análisis exploratorio

A partir del análisis exploratorio se observa que las variables con mayor relación lineal simple con pl_rade son pl_bmasse, pl_eqt, st_mass y st_rad. No obstante, al analizar las correlaciones parciales, la variable que mantiene una asociación más clara con el radio planetario es pl_bmasse, seguida de pl_eqt, sy_dist y st_mass.

También se detecta una correlación extremadamente elevada entre pl_orbper y pl_orbsmax, lo que indica un problema de multicolinealidad. Por este motivo, no se incluirán ambas variables simultáneamente en el modelo. Dado que pl_orbsmax presenta una correlación ligeramente superior con pl_rade, se conserva esta variable y se descarta pl_orbper.

Por tanto, para el ajuste inicial del modelo de regresión múltiple se considerarán como variables explicativas pl_bmasse, pl_orbsmax, pl_eqt, st_teff, st_rad, st_mass y sy_dist.

3. Ajuste del modelo lineal

fit <- lm(pl_rade ~ ., data = datos)

sfit <- summary(fit)
sfit

Call:
lm(formula = pl_rade ~ ., data = datos)

Residuals:
     Min       1Q   Median       3Q      Max
-19.0056  -2.1527  -0.6500   1.0623  15.4265

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.7440947  0.4862626   5.643 1.77e-08 ***
pl_bmasse    0.0024233  0.0001013  23.931 < 2e-16 ***
pl_orbsmax   0.0005660  0.0004823   1.174 0.240645
pl_eqt       0.0021959  0.0001427  15.384 < 2e-16 ***
st_teff     -0.0014160  0.0001614  -8.775 < 2e-16 ***
st_rad       1.1795106  0.2235132   5.277 1.38e-07 ***
st_mass      7.3293325  0.5978764  12.259 < 2e-16 ***
sy_dist     -0.0017497  0.0001327 -13.185 < 2e-16 ***

Residual standard error: 3.608 on 4358 degrees of freedom
Multiple R-squared: 0.3962, Adjusted R-squared: 0.3952
F-statistic: 408.4 on 7 and 4358 DF, p-value: < 2.2e-16

El modelo lineal múltiple se ajusta utilizando pl_rade como variable dependiente y como variables explicativas pl_bmasse, pl_orbsmax, pl_eqt, st_teff, st_rad, st_mass y sy_dist. La variable pl_orbper no se incluye debido a la multicolinealidad detectada previamente con pl_orbsmax.

El modelo global resulta estadísticamente significativo, ya que el contraste F presenta un p-valor inferior a 0.001. Por tanto, se rechaza la hipótesis nula de que todos los coeficientes del modelo sean simultáneamente iguales a cero. Esto indica que, en conjunto, las variables explicativas incluidas aportan información significativa para explicar el radio planetario.

El coeficiente de determinación múltiple es R² = 0.3962, mientras que el R² ajustado es 0.3952. Esto significa que el modelo explica aproximadamente el 39.5 % de la variabilidad observada en el radio de los exoplanetas. Aunque el modelo tiene capacidad explicativa, queda una parte importante de la variabilidad sin explicar, lo cual es razonable dada la complejidad física de los sistemas planetarios y la diversidad de tipos de exoplanetas incluidos en la muestra.

# Tabla de coeficientes
sfit$coefficients

Estimate   Std. Error    t value      Pr(>|t|)
(Intercept)   2.7440946621 0.4862626062   5.643236   1.774673e-08
pl_bmasse     0.0024233337 0.0001012619  23.931342   5.085346e-119
pl_orbsmax    0.0005660246 0.0004823243   1.173535   2.406453e-01
pl_eqt        0.0021959240 0.0001427417  15.383900   4.790326e-52
st_teff      -0.0014159586 0.0001613576  -8.775285   2.406850e-18
st_rad        1.1795105767 0.2235131926   5.277141   1.375900e-07
st_mass       7.3293325476 0.5978763670  12.258943   5.429768e-34
sy_dist      -0.0017496880 0.0001327002 -13.185274   5.889055e-39

En cuanto a los coeficientes individuales, la mayoría de variables resultan estadísticamente significativas al nivel del 5 %. En concreto, pl_bmasse, pl_eqt, st_teff, st_rad, st_mass y sy_dist presentan p-valores inferiores a 0.001, por lo que existe evidencia estadística de que están asociadas con pl_rade, manteniendo constantes el resto de variables.

La única variable no significativa es pl_orbsmax, cuyo p-valor es aproximadamente 0.241. Esto indica que, una vez controlado el efecto del resto de variables, no hay evidencia suficiente para afirmar que el semieje mayor orbital tenga un efecto lineal significativo sobre el radio planetario.

El intercepto del modelo es igual a 2.7441 radios terrestres. Este valor representa el radio planetario esperado cuando todas las variables explicativas del modelo toman simultáneamente el valor cero. Esta situación carece de interpretación física real, ya que ningún exoplaneta puede presentar simultáneamente masa nula, temperatura nula, distancia nula, etc. Por tanto, el intercepto no debe interpretarse de forma sustantiva en este contexto; forma parte de la especificación matemática del modelo.

# Intervalos de confianza para los coeficientes
confint(fit)

2.5 %       97.5 %
(Intercept)   1.7907726983  3.697416626
pl_bmasse     0.0022248089  0.002621859
pl_orbsmax   -0.0003795762  0.001511626
pl_eqt        0.0019160777  0.002475770
st_teff      -0.0017323014 -0.001099616
st_rad        0.7413110668  1.617710087
st_mass       6.1571908591  8.501474236
sy_dist      -0.0020098478 -0.001489528

coef_ci <- tidy(fit, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    significativo = ifelse(conf.low > 0 | conf.high < 0,
                           "Significativo",
                           "No significativo"),
    term = reorder(term, estimate)
  )

ggplot(coef_ci, aes(x = estimate, y = term, color = significativo)) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.8, color = "black") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.22, linewidth = 0.9) +
  geom_point(size = 3.5) +
  scale_x_continuous(
    trans = scales::pseudo_log_trans(base = 10),
    labels = scales::label_number(accuracy = 0.001)
  ) +
  scale_color_manual(values = c("Significativo" = "#E67E22", "No significativo" = "#555555")) +
  labs(
    title = "Intervalos de confianza al 95% de los coeficientes",
    x = "Estimación del coeficiente",
    y = "Variable explicativa",
    color = ""
  ) +
  theme_minimal(base_size = 13)

El coeficiente de pl_bmasse es positivo e igual a 0.00242. Esto significa que, manteniendo constantes el resto de variables, un aumento de una masa terrestre en la masa del planeta se asocia con un incremento esperado de aproximadamente 0.0024 radios terrestres en el radio planetario. Aunque el efecto unitario es pequeño, la variable es altamente significativa.

El coeficiente de pl_eqt también es positivo, con un valor aproximado de 0.00220. Esto indica que, manteniendo constantes las demás variables, un aumento de un Kelvin en la temperatura de equilibrio se asocia con un incremento esperado de 0.0022 radios terrestres en el radio planetario. En términos prácticos, un aumento de 100 K se asociaría con un incremento aproximado de 0.22 radios terrestres.

El coeficiente de st_mass es positivo y relativamente elevado, con un valor de 7.329. Esto sugiere que, manteniendo constantes el resto de variables, los planetas que orbitan estrellas de mayor masa tienden a presentar radios mayores. En concreto, un aumento de una masa solar en la masa de la estrella se asocia con un incremento esperado de aproximadamente 7.33 radios terrestres en el radio planetario.

El coeficiente de st_teff es negativo, con un valor de -0.00142. Aunque en las correlaciones simples la temperatura efectiva de la estrella presentaba una asociación positiva con pl_rade, en el modelo múltiple su efecto estimado pasa a ser negativo al controlar por el resto de variables. Esto indica que su relación con el radio planetario está condicionada por otras variables estelares, como st_mass y st_rad.

La variable sy_dist también presenta un coeficiente negativo, de aproximadamente -0.00175. Esto sugiere que, manteniendo constantes el resto de variables, los sistemas más lejanos tienden a asociarse con radios planetarios ligeramente menores en esta muestra. Esta relación puede estar influida por sesgos observacionales o por las características de los planetas mejor caracterizados.

3.1 Conclusión

Desde el punto de vista interpretativo, el modelo sugiere que el radio planetario está especialmente relacionado con la masa del planeta, la temperatura de equilibrio y las características de la estrella anfitriona. En particular, planetas más masivos, con mayor temperatura de equilibrio y orbitando estrellas de mayor masa tienden a presentar radios mayores.

No obstante, estas conclusiones deben interpretarse con cautela. El objetivo no es establecer una relación causal directa, sino identificar asociaciones estadísticas en la muestra de exoplanetas confirmados. Además, el valor de R² indica que el modelo explica una parte relevante, pero no mayoritaria, de la variabilidad del radio planetario, por lo que existen otros factores físicos no incorporados al modelo que también pueden influir en el tamaño del planeta, como la composición, la densidad, la atmósfera o la historia evolutiva del sistema.

Para identificar las variables con mayor impacto sobre el radio planetario, se ajusta también un modelo con las variables estandarizadas. Esta transformación permite comparar los coeficientes entre sí, ya que todas las variables pasan a estar medidas en la misma escala. De esta forma, el tamaño absoluto del coeficiente estandarizado refleja la importancia relativa de cada variable en el modelo.

# Modelo con variables estandarizadas
datos_estandarizados <- datos %>%
  mutate(across(everything(), scale))

fit_estandarizado <- lm(pl_rade ~ ., data = datos_estandarizados)

# Coeficientes estandarizados
coef_estandarizados <- summary(fit_estandarizado)$coefficients
coef_estandarizados

Estimate Std. Error      t value      Pr(>|t|)
(Intercept)  2.635064e-16 0.01176983 2.238830e-14  1.000000e+00
pl_bmasse    2.937203e-01 0.01227346 2.393134e+01  5.085346e-119
pl_orbsmax   1.386645e-02 0.01181596 1.173535e+00  2.406453e-01
pl_eqt       2.163914e-01 0.01406609 1.538390e+01  4.790326e-52
st_teff     -2.559607e-01 0.02916836 -8.775285e+00 2.406850e-18
st_rad       1.104647e-01 0.02093268 5.277141e+00  1.375900e-07
st_mass      4.350909e-01 0.03549171 1.225894e+01  5.429768e-34
sy_dist     -1.795005e-01 0.01361371 -1.318527e+01 5.889055e-39

coef_std_grafico <- as.data.frame(coef_estandarizados) %>%
  rownames_to_column("variable") %>%
  filter(variable != "(Intercept)") %>%
  mutate(
    significativo = ifelse(`Pr(>|t|)` < 0.05, "Significativo", "No significativo"),
    variable = reorder(variable, Estimate)
  )

ggplot(coef_std_grafico, aes(x = Estimate, y = variable, fill = significativo)) +
  geom_col() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  scale_fill_manual(values = c("Significativo" = "#E67E22", "No significativo" = "#555555")) +
  labs(
    title = "Importancia relativa de las variables",
    subtitle = "Coeficientes del modelo estandarizado",
    x = "Coeficiente estandarizado",
    y = "Variable explicativa",
    fill = ""
  ) +
  theme_minimal(base_size = 13)

# Extraemos los coeficientes sin el intercepto
coef_impacto <- coef_estandarizados[-1, , drop = FALSE]

# Los ordenamos por impacto absoluto
coef_impacto_ordenados <- coef_impacto[order(abs(coef_impacto[, "Estimate"]), decreasing = TRUE), ]

# Mostramos los 3 coeficientes con mayor impacto
coef_impacto_ordenados[1:3, ]

Estimate Std. Error   t value      Pr(>|t|)
st_mass    0.4350909 0.03549171 12.258943  5.429768e-34
pl_bmasse  0.2937203 0.01227346 23.931342  5.085346e-119
st_teff   -0.2559607 0.02916836 -8.775285  2.406850e-18

Según los coeficientes estandarizados, las tres variables con mayor impacto relativo sobre pl_rade son st_mass, pl_bmasse y st_teff. La variable con mayor efecto es st_mass, con un coeficiente positivo de 0.435, lo que indica que los planetas que orbitan estrellas de mayor masa tienden a presentar radios mayores, manteniendo constantes el resto de variables.

En segundo lugar aparece pl_bmasse, también con efecto positivo, por lo que los planetas más masivos tienden a tener mayor radio. Por último, st_teff presenta un coeficiente negativo, lo que sugiere que, al controlar por el resto de variables del modelo, una mayor temperatura efectiva de la estrella se asocia con un menor radio planetario. Las tres variables son estadísticamente significativas, ya que sus p-valores son inferiores a 0.001.

3.2 Recomendaciones

A partir de los resultados del modelo, se formulan recomendaciones orientadas a identificar qué tipo de sistemas planetarios presentan, en promedio, mayores radios planetarios. Estas recomendaciones se basan en el signo, la magnitud, la precisión y la significatividad estadística de los coeficientes estimados, manteniendo siempre constantes el resto de variables del modelo.

Priorizar sistemas con estrellas de mayor masa (st_mass). Esta es la variable con mayor impacto relativo según los coeficientes estandarizados (0.435). Su coeficiente es positivo, significativo (p < 0.001) e interpretable: un aumento de una masa solar se asocia con un incremento esperado de aproximadamente 7.33 radios terrestres en el radio planetario.

Seleccionar planetas con mayor masa (pl_bmasse). Con un coeficiente estandarizado de 0.232 y un coeficiente sin estandarizar de 0.00242 (p < 0.001), la masa planetaria presenta una relación positiva y significativa con el radio. Aunque el efecto unitario es pequeño, la amplitud reducida del intervalo de confianza indica una estimación precisa.

Considerar planetas con mayor temperatura de equilibrio (pl_eqt). Su coeficiente es positivo (0.00220) y estadísticamente significativo (p < 0.001). Un aumento de 100 K en la temperatura de equilibrio se asocia con aproximadamente 0.22 radios terrestres adicionales.

Cautela con st_teff, sy_dist y pl_orbsmax. La temperatura efectiva de la estrella presenta un coeficiente negativo en el modelo múltiple (-0.00142), la distancia al sistema también muestra un coeficiente negativo y significativo, y pl_orbsmax no resulta estadísticamente significativa (p ≈ 0.241).

En conjunto, el modelo sugiere que los sistemas con estrellas masivas y planetas igualmente masivos y calientes son los mejores candidatos para encontrar exoplanetas de mayor tamaño. No obstante, estas recomendaciones deben interpretarse como asociaciones estadísticas y no como relaciones causales, dado que el modelo explica únicamente el 39.5 % de la variabilidad del radio planetario.

4. Bondad del ajuste

# R^2
sfit$r.squared

# R^2 ajustado
sfit$adj.r.squared

# CV-sigma
cv_sigma <- sigma(fit) / mean(datos$pl_rade) * 100
cv_sigma

[1] 0.3961529
[1] 0.395183
[1] 80.63467

El modelo presenta un R² de 0.3962 y un R² ajustado de 0.3952, por lo que explica aproximadamente el 39.5 % de la variabilidad del radio planetario. El contraste global del modelo es significativo, ya que el p-valor asociado al estadístico F es inferior a 0.001, lo que indica que las variables explicativas aportan información conjunta sobre pl_rade.

Sin embargo, el valor de CV_sigma es elevado, aproximadamente 80.63 %. Este indicador compara el error residual del modelo con la media de la variable dependiente. Un valor tan alto indica que el error residual es considerable en relación con el tamaño medio de pl_rade. Por tanto, aunque el modelo es estadísticamente significativo, su capacidad predictiva es moderada y podría mejorarse mediante transformaciones o modelos alternativos.

5. Análisis de los residuos

Una vez ajustado el modelo, no basta con mirar si los coeficientes son significativos o si el R² parece aceptable. En regresión lineal hay que comprobar si el modelo cumple sus supuestos básicos. Para eso se analizan los residuos, que son la diferencia entre el valor observado y el valor predicho por el modelo:

residuo_i = y_i - \hat{y}_i

Si el modelo está bien especificado, los residuos deberían comportarse como un “ruido” alrededor de cero: sin patrones claros, con varianza aproximadamente constante, con una distribución razonablemente normal y sin observaciones que dominen el ajuste. Justo por eso R genera cuatro gráficos de diagnóstico cuando ejecutamos plot(fit).

plot(fit)

5.1 ¿Qué representa cada gráfico de diagnóstico?

1) Residuals vs Fitted

Este gráfico coloca en el eje X los valores ajustados por el modelo, es decir, las predicciones \hat{y}, y en el eje Y los residuos. Sirve para comprobar principalmente dos cosas: si la relación entre las variables puede aproximarse razonablemente con una forma lineal y si la dispersión de los errores se mantiene más o menos constante.

Lo ideal sería ver una nube de puntos sin forma clara, centrada alrededor de la línea horizontal de cero. En este caso no ocurre eso: aparece un patrón curvado y una forma de abanico. Esto sugiere que el modelo lineal con las variables originales no está capturando bien toda la estructura de los datos. Dicho de otra forma, hay señales de falta de linealidad y de posible heterocedasticidad, porque la variabilidad de los residuos cambia según el nivel de predicción.

2) Q-Q Residuals

El gráfico Q-Q compara los cuantiles de los residuos estandarizados con los cuantiles que esperaríamos si siguieran una distribución normal. Es una forma visual de evaluar el supuesto de normalidad de los residuos.

Si los residuos fueran aproximadamente normales, los puntos caerían cerca de la diagonal. En el gráfico se observa una desviación clara, sobre todo en las colas: los puntos de los extremos se separan bastante de la línea. Esto indica que los residuos tienen valores extremos y que la normalidad no se cumple bien. Este resultado es coherente con lo que ya se veía desde el histograma de pl_rade: la variable dependiente original era asimétrica y contenía valores grandes.

3) Scale-Location

El gráfico Scale-Location representa la raíz de los residuos estandarizados absolutos frente a los valores ajustados. Aunque suene más técnico, su objetivo es bastante directo: comprobar si la varianza de los residuos es constante.

Lo deseable sería una banda de puntos con una dispersión parecida a lo largo de todo el eje X y una línea roja más o menos horizontal. En este caso la línea roja cambia claramente y la dispersión aumenta en determinadas zonas. Esto refuerza la idea de que hay heterocedasticidad: el modelo no comete errores con la misma variabilidad para todos los tamaños de exoplanetas.

4) Residuals vs Leverage

Este gráfico sirve para localizar observaciones que pueden tener mucha influencia sobre el modelo. Aquí aparecen dos conceptos: leverage e influencia. El leverage mide si una observación tiene valores raros o extremos en las variables explicativas; la influencia mide si, además de ser rara, cambia de forma importante el ajuste del modelo.

Las curvas discontinuas corresponden a referencias de distancia de Cook. Si una observación aparece con leverage alto y residuo grande, puede tener un efecto importante sobre los coeficientes. En el gráfico se ven algunas observaciones destacadas, por lo que tiene sentido completar el análisis con la distancia de Cook de forma explícita.

En conjunto, estos cuatro gráficos indican que el modelo lineal inicial no cumple completamente los supuestos clásicos. Esto no significa que el modelo no sirva para explorar asociaciones, pero sí obliga a interpretar los p-valores, los intervalos de confianza y la capacidad predictiva con cautela.

5.2 Histograma de residuos

# Extraemos residuos
 e <- residuals(fit)

# Histograma de residuos
hist(e)

Aquí sí corresponde mostrar una imagen, porque hist(e) genera un gráfico. El histograma permite ver cómo se distribuyen los errores del modelo. La mayor parte de los residuos se concentran cerca de cero, lo cual es positivo: significa que muchas predicciones no se alejan demasiado del valor observado. Sin embargo, la distribución no es perfectamente simétrica y aparecen colas con valores extremos, especialmente hacia residuos negativos grandes y positivos altos.

Esto es importante porque en el modelo lineal clásico se asume que los errores tienen una distribución aproximadamente normal. Si las colas son muy pesadas o hay asimetría, los contrastes de significación y los intervalos de confianza pueden perder fiabilidad. Por eso el histograma no se interpreta de forma aislada: se complementa con el Q-Q plot y con el test de Shapiro-Wilk.

5.3 Contraste de normalidad: Shapiro-Wilk

# Contraste de normalidad de Shapiro-Wilk
shapiro.test(e)

Shapiro-Wilk normality test

data:  e
W = 0.93023, p-value < 2.2e-16

El test de Shapiro-Wilk contrasta la hipótesis nula de que los residuos proceden de una distribución normal. En este caso el p-valor es inferior a 0.001, por lo que se rechaza la normalidad. Este resultado confirma formalmente lo que ya se veía en el Q-Q plot y en el histograma: los residuos no se ajustan bien a una distribución normal.

En muestras grandes, como ocurre aquí, los tests de normalidad son muy sensibles y pueden rechazar la normalidad incluso ante desviaciones pequeñas. Aun así, en este proyecto la evidencia visual también muestra desviaciones claras, de modo que no conviene ignorar el resultado.

5.4 Contraste de homocedasticidad: Breusch-Pagan

# Contraste de homocedasticidad de Breusch-Pagan
bptest(fit)

studentized Breusch-Pagan test

data:  fit
BP = 1369.1, df = 7, p-value < 2.2e-16

La homocedasticidad significa que la varianza de los errores debería ser aproximadamente constante para todos los niveles de predicción. En términos sencillos: el modelo no debería equivocarse mucho más para unos planetas que para otros solo porque el radio predicho sea mayor o menor.

El test de Breusch-Pagan tiene como hipótesis nula que la varianza de los residuos es constante. Como el p-valor es inferior a 0.001, se rechaza esa hipótesis. Por tanto, hay evidencia de heterocedasticidad. Esto encaja con los gráficos Residuals vs Fitted y Scale-Location, donde la dispersión de los residuos cambia según los valores ajustados.

5.5 Contraste de independencia: Durbin-Watson

# Contraste de incorrelación de Durbin-Watson
dwtest(fit, alternative = "two.sided")

Durbin-Watson test

data:  fit
DW = 1.107, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0

El test de Durbin-Watson se utiliza para comprobar si existe autocorrelación en los residuos. La idea es ver si los errores del modelo están relacionados entre sí en lugar de comportarse como errores independientes.

El estadístico Durbin-Watson suele interpretarse tomando el valor 2 como referencia: valores cercanos a 2 indican ausencia de autocorrelación; valores claramente inferiores a 2 sugieren autocorrelación positiva; valores superiores a 2 sugieren autocorrelación negativa. Aquí se obtiene DW = 1.107 con p-valor inferior a 0.001, por lo que se rechaza la hipótesis de ausencia de autocorrelación. Esto indica que los residuos presentan autocorrelación positiva.

5.6 Valores sobreinfluyentes y distancia de Cook

Una observación sobreinfluyente es una observación que, si se eliminara del conjunto de datos, podría cambiar de forma apreciable los coeficientes estimados del modelo. No es simplemente un dato raro: es un dato con capacidad real de mover el ajuste.

Para medir esto se utiliza la distancia de Cook. Esta medida combina información sobre el residuo y sobre el leverage de una observación. Una observación puede tener un residuo grande, pero si no cambia el ajuste global quizá no sea tan preocupante. También puede tener leverage alto, pero si está bien alineada con el modelo quizá tampoco distorsione demasiado. La distancia de Cook detecta los casos en los que la observación es relevante para el ajuste completo.

# Distancia de Cook
cook <- cooks.distance(fit)
n <- nrow(datos)
umbral <- 4 / n

# Número de observaciones sobreinfluyentes
sum(cook > umbral)

# Proporción sobre el total
mean(cook > umbral) * 100

[1] 265
[1] 6.069629

# Gráfico de distancias de Cook
plot(cook, type = "h", main = "Distancia de Cook por observación",
     ylab = "Distancia de Cook", xlab = "Índice")
abline(h = umbral, col = "red", lty = 2)

El umbral utilizado es 4/n, una regla práctica habitual para señalar observaciones potencialmente influyentes. En este caso se detectan 265 observaciones sobreinfluyentes, lo que representa aproximadamente el 6.1 % de la muestra. No es la mayoría de los datos, pero sí es una cantidad suficiente como para tenerla en cuenta.

La conclusión no debería ser eliminar automáticamente esas observaciones. En un contexto astronómico, algunos valores extremos pueden ser reales y físicamente interesantes. Lo correcto es interpretarlas como una advertencia: el ajuste lineal está condicionado por una parte de la muestra y conviene probar transformaciones o modelos más robustos antes de extraer conclusiones demasiado fuertes.

Por tanto, el análisis de residuos deja una conclusión clara: el modelo lineal con variables originales es útil como primera aproximación, pero no cumple plenamente normalidad, homocedasticidad, independencia ni ausencia de observaciones influyentes. Esta es la razón estadística que justifica la ampliación posterior con transformación logarítmica.

6. Selección de variables

# Selección automática de variables mediante AIC
fit_step <- step(fit, direction = "backward")

# Resumen del modelo seleccionado
sfit_step <- summary(fit_step)
sfit_step

Start: AIC=11213.39
pl_rade ~ pl_bmasse + pl_orbsmax + pl_eqt + st_teff + st_rad + st_mass + sy_dist

             Df Sum of Sq   RSS   AIC
- pl_orbsmax  1      17.9 56761 11213
<none>                    56743 11213
- st_rad      1     362.6 57105 11239
- st_teff     1    1002.6 57745 11288
- st_mass     1    1956.7 58699 11359
- sy_dist     1    2263.6 59006 11382
- pl_eqt      1    3081.5 59824 11442
- pl_bmasse   1    7456.9 64200 11750

Step: AIC=11212.77
pl_rade ~ pl_bmasse + pl_eqt + st_teff + st_rad + st_mass + sy_dist

Residual standard error: 3.609 on 4359 degrees of freedom
Multiple R-squared: 0.396, Adjusted R-squared: 0.3951
F-statistic: 476.2 on 6 and 4359 DF, p-value: < 2.2e-16

Se aplica el procedimiento de selección hacia atrás basado en el criterio AIC. Este algoritmo parte del modelo completo y va eliminando variables cuando su eliminación mejora el AIC. El objetivo es obtener un modelo más simple, manteniendo una capacidad explicativa adecuada.

# R2 del modelo simplificado
sfit_step$r.squared

# R2 ajustado del modelo simplificado
sfit_step$adj.r.squared

# AIC del modelo completo y del modelo simplificado
AIC(fit, fit_step)

# CV-sigma del modelo simplificado
cv_sigma_step <- sigma(fit_step) / mean(datos$pl_rade) * 100
cv_sigma_step

[1] 0.3959621
[1] 0.3951306

         df      AIC
fit       9 23605.56
fit_step  8 23604.94

[1] 80.63816

comparacion_modelos <- data.frame(
  Modelo = c("Modelo general", "Modelo AIC"),
  R2 = c(sfit$r.squared, sfit_step$r.squared),
  R2_ajustado = c(sfit$adj.r.squared, sfit_step$adj.r.squared),
  AIC = c(AIC(fit), AIC(fit_step)),
  CV_sigma = c(
    sigma(fit) / mean(datos$pl_rade) * 100,
    sigma(fit_step) / mean(datos$pl_rade) * 100
  )
)

comparacion_modelos

Modelo        R2 R2_ajustado      AIC CV_sigma
1 Modelo general 0.3961529   0.3951830 23605.56 80.63467
2    Modelo AIC 0.3959621   0.3951306 23604.94 80.63816

Ambos modelos presentan métricas de bondad del ajuste prácticamente idénticas. El modelo AIC elimina pl_orbsmax, que no resultaba significativa, y obtiene un AIC ligeramente inferior (23604.94 frente a 23605.56), por lo que se prefiere el modelo simplificado por ser más parsimonioso sin perder capacidad explicativa.

# Gráficos de diagnóstico del modelo AIC
plot(fit_step)

# Residuos del modelo AIC
e_step <- residuals(fit_step)

# Histograma
hist(e_step)

# Normalidad
shapiro.test(e_step)

# Homocedasticidad
bptest(fit_step)

# Autocorrelación
dwtest(fit_step, alternative = "two.sided")

Shapiro-Wilk: W = 0.93015, p-value < 2.2e-16
Breusch-Pagan: BP = 1371.9, df = 6, p-value < 2.2e-16
Durbin-Watson: DW = 1.1061, p-value < 2.2e-16

test_residuos <- data.frame(
  Modelo = c("Modelo general", "Modelo AIC"),
  Shapiro_pvalor = c(
    shapiro.test(residuals(fit))$p.value,
    shapiro.test(residuals(fit_step))$p.value
  ),
  Breusch_Pagan_pvalor = c(
    bptest(fit)$p.value,
    bptest(fit_step)$p.value
  ),
  Durbin_Watson_pvalor = c(
    dwtest(fit, alternative = "two.sided")$p.value,
    dwtest(fit_step, alternative = "two.sided")$p.value
  )
)

test_residuos

Modelo Shapiro_pvalor Breusch_Pagan_pvalor Durbin_Watson_pvalor
1 Modelo general   3.988008e-41         1.874945e-291        9.399645e-192
2    Modelo AIC   3.836369e-41         2.998720e-293        3.819913e-192

Sin embargo, ninguno de los dos modelos permite realizar inferencia estadística de forma completamente válida, ya que ambos incumplen los supuestos básicos de normalidad, homocedasticidad e independencia de los residuos. Los p-valores y los intervalos de confianza deben interpretarse con cautela, ya que su validez depende del cumplimiento de dichos supuestos.

7. Conclusiones finales

En este trabajo se ha analizado un conjunto de datos de exoplanetas confirmados con el objetivo de estudiar qué variables físicas y orbitales se relacionan con el radio planetario (pl_rade). Para ello, se realizó primero un análisis exploratorio de los datos, posteriormente se ajustó un modelo de regresión lineal múltiple y, finalmente, se evaluó la calidad del ajuste y el cumplimiento de los supuestos del modelo.

El análisis exploratorio mostró que la variable dependiente presenta una distribución asimétrica y con valores extremos, algo esperable en una base de datos astronómica formada por planetas de características muy diferentes. Además, se detectó una relación muy elevada entre pl_orbper y pl_orbsmax, lo que indicaba un problema de multicolinealidad. Por este motivo, se decidió conservar pl_orbsmax y excluir pl_orbper del modelo inicial.

El modelo lineal múltiple resultó globalmente significativo, lo que indica que las variables explicativas consideradas aportan información relevante para explicar el radio planetario. En particular, variables como pl_bmasse, pl_eqt, st_teff, st_rad, st_mass y sy_dist resultaron estadísticamente significativas. La única variable que no mostró un efecto significativo en el modelo fue pl_orbsmax.

Desde el punto de vista interpretativo, los resultados indican que el radio planetario está especialmente relacionado con la masa del planeta y con características de la estrella anfitriona. El modelo con variables estandarizadas permitió comparar la importancia relativa de los predictores, destacando principalmente st_mass, pl_bmasse y st_teff como las variables con mayor impacto relativo sobre pl_rade.

No obstante, la capacidad explicativa del modelo es moderada. El valor de R² ajustado indica que el modelo explica aproximadamente el 39.5 % de la variabilidad del radio planetario, por lo que una parte importante de la variación queda sin explicar. Esto es razonable, ya que el tamaño de un exoplaneta depende de muchos factores físicos que no están incluidos en el modelo, como la composición interna, la densidad, la atmósfera, la edad del sistema o los procesos de formación y evolución planetaria.

Además, el análisis de residuos mostró que el modelo no cumple completamente los supuestos clásicos de regresión lineal. En particular, se detectaron problemas de normalidad, heterocedasticidad y autocorrelación en los residuos. Por tanto, aunque el modelo es útil como primera aproximación explicativa, sus resultados deben interpretarse con cautela y no deben considerarse como una herramienta predictiva definitiva.

La selección automática de variables mediante AIC permitió obtener un modelo más simple, manteniendo una capacidad explicativa similar. Sin embargo, los problemas detectados en los residuos indican que la mejora del modelo no debería centrarse únicamente en eliminar variables, sino también en considerar transformaciones o enfoques alternativos.

8. Ampliación: modelo con transformación logarítmica

Como ampliación del análisis, se plantea un modelo alternativo basado en transformaciones logarítmicas. Esta propuesta se justifica por la elevada asimetría observada en algunas variables del conjunto de datos, especialmente en pl_rade, pl_bmasse, pl_orbsmax y sy_dist. Estas variables presentan valores muy dispersos y algunos valores extremos, lo que puede afectar al cumplimiento de los supuestos del modelo lineal clásico.

La transformación logarítmica puede ayudar a reducir la influencia de valores extremos, estabilizar parcialmente la varianza y mejorar la relación lineal entre la variable dependiente y algunas variables explicativas. Por este motivo, se ajusta un modelo donde se transforma logarítmicamente el radio planetario y algunas variables estrictamente positivas.

# Creamos una nueva base con transformaciones logarítmicas
datos_log <- datos %>%
  mutate(
    log_pl_rade = log(pl_rade),
    log_pl_bmasse = log(pl_bmasse),
    log_pl_orbsmax = log(pl_orbsmax),
    log_sy_dist = log(sy_dist)
  )

# Ajuste del modelo logarítmico
fit_log <- lm(
  log_pl_rade ~ log_pl_bmasse + log_pl_orbsmax + pl_eqt +
    st_teff + st_rad + st_mass + log_sy_dist,
  data = datos_log
)

# Resumen del modelo
summary(fit_log)

Call:
lm(formula = log_pl_rade ~ log_pl_bmasse + log_pl_orbsmax + pl_eqt +
    st_teff + st_rad + st_mass + log_sy_dist, data = datos_log)

Residuals:
     Min       1Q   Median       3Q      Max
-3.08251 -0.12595  0.00827  0.13445  1.52787

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)      3.116e-02  6.168e-02   0.505 0.613458
log_pl_bmasse    3.772e-01  3.083e-03 122.354 < 2e-16 ***
log_pl_orbsmax  -5.114e-03  1.024e-02  -0.499 0.617532
pl_eqt          -2.154e-05  2.475e-05  -0.870 0.384298
st_teff         -2.454e-05  1.655e-05  -1.483 0.138181
st_rad           7.574e-02  2.182e-02   3.470 0.000525 ***
st_mass          2.322e-01  5.762e-02   4.029 5.69e-05 ***
log_sy_dist     -5.867e-03  5.396e-03  -1.087 0.276976

Residual standard error: 0.344 on 4358 degrees of freedom
Multiple R-squared: 0.8347, Adjusted R-squared: 0.8345
F-statistic: 3144 on 7 and 4358 DF, p-value: < 2.2e-16

En este modelo, la variable dependiente es log_pl_rade, es decir, el logaritmo del radio planetario. Por tanto, la interpretación de los coeficientes cambia respecto al modelo lineal original. En el caso de las variables también transformadas en logaritmos, como log_pl_bmasse, log_pl_orbsmax y log_sy_dist, los coeficientes pueden interpretarse aproximadamente como elasticidades.

8.1 Comparación de la bondad del ajuste

comparacion_modelos <- tibble(
  Modelo = c("Modelo general", "Modelo AIC", "Modelo logarítmico"),
  R2 = c(summary(fit)$r.squared, summary(fit_step)$r.squared, summary(fit_log)$r.squared),
  R2_ajustado = c(summary(fit)$adj.r.squared, summary(fit_step)$adj.r.squared, summary(fit_log)$adj.r.squared),
  AIC = c(AIC(fit), AIC(fit_step), AIC(fit_log)),
  Sigma_residual = c(sigma(fit), sigma(fit_step), sigma(fit_log))
)

knitr::kable(comparacion_modelos, digits = 4,
             caption = "Comparación de la bondad del ajuste entre modelos")

Modelo              R2     R2_ajustado      AIC       Sigma_residual
Modelo general      0.3962 0.3952       23605.559     3.6084
Modelo AIC          0.3960 0.3951       23604.939     3.6085
Modelo logarítmico  0.8347 0.8345        3081.181     0.3440

El modelo logarítmico se compara con los modelos anteriores para comprobar si la transformación mejora la capacidad explicativa, reduce el AIC o disminuye la dispersión de los residuos. En particular, un mayor valor de R² ajustado, un menor AIC y una menor desviación típica residual indicarían una mejora relativa del ajuste.

No obstante, debe tenerse en cuenta que el modelo logarítmico no está expresado en la misma escala que el modelo original, ya que la variable dependiente pasa de ser pl_rade a log_pl_rade. Por ello, la comparación debe interpretarse con cautela, especialmente en el caso de la desviación típica residual.

8.2 Intervalos de confianza del modelo logarítmico

# Intervalos de confianza del modelo logarítmico
confint(fit_log)

2.5 %        97.5 %
(Intercept)     -8.975845e-02  1.520730e-01
log_pl_bmasse    3.711205e-01  3.832072e-01
log_pl_orbsmax  -2.518869e-02  1.496151e-02
pl_eqt          -7.007009e-05  2.699273e-05
st_teff         -5.698024e-05  7.903804e-06
st_rad           3.294918e-02  1.185225e-01
st_mass          1.192021e-01  3.451204e-01
log_sy_dist     -1.644491e-02  4.711703e-03

A partir de los intervalos de confianza al 95 %, se observa que log_pl_bmasse, st_rad y st_mass presentan efectos estadísticamente significativos, ya que sus intervalos no contienen el valor 0. En concreto, log_pl_bmasse muestra un intervalo claramente positivo, lo que confirma que la masa planetaria sigue siendo la variable más relevante para explicar el radio planetario en el modelo logarítmico.

En cambio, los intervalos de log_pl_orbsmax, pl_eqt, st_teff y log_sy_dist contienen el valor 0, por lo que no se puede afirmar que estas variables tengan un efecto estadísticamente significativo sobre el logaritmo del radio planetario al nivel de confianza del 95 %.

8.3 Diagnóstico de residuos del modelo logarítmico

# Gráficos diagnósticos del modelo logarítmico
par(mfrow = c(2, 2))
plot(fit_log)
par(mfrow = c(1, 1))

Los gráficos diagnósticos permiten comprobar visualmente si la transformación logarítmica mejora el comportamiento de los residuos. En concreto, se analiza si los residuos presentan una distribución más cercana a la normalidad, si la varianza es más constante y si se reduce la influencia de observaciones extremas.

# Contraste de normalidad de los residuos
shapiro.test(residuals(fit_log))

# Contraste de homocedasticidad
bptest(fit_log)

# Contraste de autocorrelación
dwtest(fit_log)

# Valores influyentes mediante distancia de Cook
cook_log <- cooks.distance(fit_log)
umbral_cook_log <- 4 / nrow(datos_log)

influyentes_log <- which(cook_log > umbral_cook_log)

length(influyentes_log)
head(influyentes_log)

Shapiro-Wilk: W = 0.8486, p-value < 2.2e-16
Breusch-Pagan: BP = 499.68, df = 7, p-value < 2.2e-16
Durbin-Watson: DW = 1.4986, p-value < 2.2e-16

[1] 257
[1]  3  4  6 13 19 24

Los contrastes estadísticos muestran que el modelo logarítmico mejora el ajuste, pero sigue presentando problemas en los supuestos básicos. El test de Shapiro-Wilk rechaza la normalidad de los residuos, el test de Breusch-Pagan rechaza la hipótesis de homocedasticidad y el test de Durbin-Watson detecta autocorrelación positiva en los residuos.

Por último, el análisis mediante distancia de Cook identifica 257 observaciones potencialmente influyentes. Esto confirma que, aunque la transformación logarítmica reduce parte de la asimetría y mejora notablemente la bondad del ajuste, todavía existen observaciones extremas o influyentes que afectan al comportamiento del modelo.

8.4 Conclusión de la ampliación

La transformación logarítmica mejora de forma clara la capacidad explicativa del modelo. Mientras que el modelo general y el modelo seleccionado mediante AIC explicaban aproximadamente el 39.5 % de la variabilidad del radio planetario, el modelo logarítmico alcanza un R² ajustado de 0.8345. Además, el AIC se reduce notablemente, pasando de valores cercanos a 23605 en los modelos anteriores a 3081.18 en el modelo logarítmico.

También se reduce considerablemente la desviación típica residual, que pasa de aproximadamente 3.61 en los modelos con la variable original a 0.344 en el modelo logarítmico. No obstante, esta comparación debe interpretarse con cautela, ya que el modelo logarítmico no predice directamente el radio planetario en radios terrestres, sino su logaritmo.

En cuanto a los coeficientes del modelo logarítmico, las variables estadísticamente significativas son log_pl_bmasse, st_rad y st_mass. La variable con mayor importancia es log_pl_bmasse, cuyo coeficiente positivo indica que los planetas con mayor masa tienden a presentar un mayor radio, manteniendo constantes el resto de variables.

Sin embargo, el diagnóstico de residuos muestra que la transformación logarítmica no resuelve completamente las limitaciones del modelo. Por tanto, el modelo logarítmico puede considerarse una mejora importante respecto al modelo lineal original en términos de capacidad explicativa y ajuste global, pero no cumple plenamente los supuestos clásicos de regresión lineal. En conjunto, la transformación logarítmica permite capturar mejor la relación entre las variables astronómicas y el radio planetario, pero la complejidad y heterogeneidad de los datos sugieren que podrían ser necesarios modelos más flexibles, como regresión robusta, modelos no lineales o modelos específicos para poblaciones diferenciadas de exoplanetas.