Análisis de Masa y Luminosidad Estelar

Este análisis utiliza datos de temperatura efectiva (\(T_{\text{eff}}\)) para calcular propiedades estelares como masa y luminosidad. Además, se aplica la transformación de Box-Cox para mejorar las relaciones lineales y se realizan predicciones sobre nuevos datos.

Carga de Librerías y Datos

Se cargan las librerías necesarias y los datos provenientes de un archivo .ods. Este conjunto de datos incluye información sobre la temperatura efectiva (\(T_{\text{eff}}\)) de diferentes estrellas.

library(ggplot2)
library(MASS)
library(readODS)

file_path <- "/home/cesar/Documentos/Masas.ods" # Cambiar la ruta acorde a la ubicación en tu computadora
data <- read_ods(file_path, sheet = 1)
head(data)

## # A tibble: 6 × 10
##   bestobjid specobjid plate   mjd fiberid T_eff   Fe_H Log_g   snr flag 
##       <dbl>     <dbl> <dbl> <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl> <chr>
## 1   1.24e18   3.28e17   291 51928     321 5607. -0.595  4.22  40.2 nnnnn
## 2   1.24e18   3.28e17   291 51928     374 5551. -0.277  4.10  35.8 nnBnn
## 3   1.24e18   3.28e17   291 51928     426 5802. -0.785  4.20  34.5 nnnnn
## 4   1.24e18   3.28e17   291 51928     471 6206. -1.48   4.11  52.5 nnnnn
## 5   1.24e18   3.64e17   323 51615       2 4235. -1.23   4.30  30.9 nnngn
## 6   1.24e18   3.64e17   323 51615       5 4177. -1.50   4.48  30.4 nnngn

Cálculo de Masa Estelar

La masa estelar se calcula en función de la temperatura efectiva usando una relación proporcional: - Si \(T_{\text{eff}} < T_{\odot}\): \(M \propto \left(\frac{T_{\text{eff}}}{T_{\odot}}\right)^2\) - Si \(T_{\text{eff}} \geq T_{\odot}\): \(M \propto \left(\frac{T_{\text{eff}}}{T_{\odot}}\right)^4\)

T_sun <- 5778  # Temperatura en Kelvin

# Función para calcular masa estelar
calculate_mass <- function(T_eff) {
  if (T_eff < T_sun) {
    return((T_eff / T_sun)^2)
  } else {
    return((T_eff / T_sun)^4)
  }
}

# Aplicar función a los datos
data$Masa_estelar <- sapply(data$T_eff, calculate_mass)
head(data[c("T_eff", "Masa_estelar")])

## # A tibble: 6 × 2
##   T_eff Masa_estelar
##   <dbl>        <dbl>
## 1 5607.        0.942
## 2 5551.        0.923
## 3 5802.        1.02 
## 4 6206.        1.33 
## 5 4235.        0.537
## 6 4177.        0.523

Cálculo de Luminosidad Estelar

La luminosidad estelar (\(L\)) se calcula en función de la masa (\(M\)) y la temperatura efectiva (\(T_{\text{eff}}\)). La relación utilizada considera que:

El radio estelar (\(R\)) es proporcional a \(M^{0.8}\).
La luminosidad es proporcional a \(R^2 \cdot \left(\frac{T_{\text{eff}}}{T_{\odot}}\right)^4\).

Por lo tanto, la fórmula es: \[ L \propto R^2 \cdot \left(\frac{T_{\text{eff}}}{T_{\odot}}\right)^4, \quad \text{donde} \quad R \propto M^{0.8}. \]

# Función para calcular luminosidad
calculate_luminosity <- function(mass, T_eff) {
  R <- mass^0.8 
  L <- R^2 * (T_eff / T_sun)^4
  return(L)
}

# Aplicar función a los datos
data$Luminosidad <- mapply(calculate_luminosity, data$Masa_estelar, data$T_eff)
head(data[c("T_eff", "Masa_estelar", "Luminosidad")])

## # A tibble: 6 × 3
##   T_eff Masa_estelar Luminosidad
##   <dbl>        <dbl>       <dbl>
## 1 5607.        0.942      0.806 
## 2 5551.        0.923      0.750 
## 3 5802.        1.02       1.04  
## 4 6206.        1.33       2.10  
## 5 4235.        0.537      0.107 
## 6 4177.        0.523      0.0968

# Relación Masa - T_eff
ggplot(data, aes(x = T_eff, y = Masa_estelar)) +
  geom_point(color = "blue", alpha = 0.7) +
  labs(title = "Relación Masa - T_eff", x = "Temperatura Efectiva (K)", y = "Masa Estelar") +
  theme_minimal()

# Relación Luminosidad - T_eff
ggplot(data, aes(x = T_eff, y = Luminosidad)) +
  geom_point(color = "red", alpha = 0.7) +
  labs(title = "Relación Luminosidad - T_eff", x = "Temperatura Efectiva (K)", y = "Luminosidad") +
  theme_minimal()

# Relación entre Masa Estelar y Luminosidad
ggplot(data, aes(x = Masa_estelar, y = Luminosidad)) +
  geom_point(color = "purple", alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", color = "orange", linetype = "dashed", se = FALSE) + # Línea de tendencia
  labs(
    title = "Relación entre Masa Estelar y Luminosidad",
    x = "Masa Estelar (M/M_⊙)",
    y = "Luminosidad (L/L_⊙)"
  ) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

# M_⊙ := Masa del Sol
# L_⊙ := Luminosidad del Sol

# Histograma y densidad de la masa estelar
ggplot(data, aes(x = Masa_estelar)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.1, fill = "blue", color = "black", alpha = 0.7) +
  labs(
    title = "Distribución de la Masa Estelar",
    x = "Masa Estelar (M/M_sun)",
    y = "Densidad"
  ) +
  theme_minimal()

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Histograma y densidad de la luminosidad
ggplot(data, aes(x = Luminosidad)) +
  geom_histogram(aes(y = ..density..), binwidth = 1.5, fill = "green", color = "black", alpha = 0.7)+
  labs(
    title = "Distribución de la Luminosidad",
    x = "Luminosidad (L/L_sun)",
    y = "Densidad"
  ) +
  theme_minimal()

## Transformaciones Box-Cox

Se aplica la transformación de Box-Cox para mejorar la linealidad de las relaciones entre variables.

# Función para calcular la verosimilitud logarítmica de Box-Cox
log_likelihood_boxcox <- function(y, lambda) {
  n <- length(y)
  if (lambda == 0) {
    # Caso especial: logaritmo
    y_trans <- log(y)
  } else {
    # Transformación general
    y_trans <- (y^lambda - 1) / lambda
  }
  log_likelihood <- -n / 2 * log(var(y_trans)) + (lambda - 1) * sum(log(y))
  return(log_likelihood)
}

# Encontrar el mejor lambda
find_optimal_lambda <- function(y, lambda_seq) {
  log_likelihoods <- sapply(lambda_seq, log_likelihood_boxcox, y = y)
  optimal_lambda <- lambda_seq[which.max(log_likelihoods)]
  return(list(optimal_lambda = optimal_lambda, log_likelihoods = log_likelihoods))
}

# Aplicar la transformación Box-Cox
apply_boxcox <- function(y, lambda) {
  if (lambda == 0) {
    return(log(y))
  } else {
    return((y^lambda - 1) / lambda)
  }
}

# Secuencia de valores lambda usando método de Grid Search
lambda_seq <- seq(-2, 2, 0.1)

# Aplicar para Masa Estelar
result_masa <- find_optimal_lambda(data$Masa_estelar, lambda_seq)
lambda_masa <- result_masa$optimal_lambda
cat("Lambda óptimo para Masa Estelar:", lambda_masa, "\n")

## Lambda óptimo para Masa Estelar: -0.8

data$Masa_estelar_BoxCox <- apply_boxcox(data$Masa_estelar, lambda_masa)

# Aplicar para Luminosidad
result_luminosidad <- find_optimal_lambda(data$Luminosidad, lambda_seq)
lambda_luminosidad <- result_luminosidad$optimal_lambda
cat("Lambda óptimo para Luminosidad:", lambda_luminosidad, "\n")

## Lambda óptimo para Luminosidad: -0.2

data$Luminosidad_BoxCox <- apply_boxcox(data$Luminosidad, lambda_luminosidad)

# Función general para calcular log-verosimilitud e intervalo de confianza
calculate_ci <- function(log_likelihoods, lambda_seq) {
  max_ll <- max(log_likelihoods)
  cutoff <- max_ll - 0.5 * qchisq(0.95, df = 1)
  ci_indices <- which(log_likelihoods >= cutoff)
  ci_lambda <- range(lambda_seq[ci_indices])
  return(list(cutoff = cutoff, ci_lambda = ci_lambda))
}

# Función para graficar log-verosimilitud
plot_log_likelihood <- function(log_likelihoods, lambda_seq, optimal_lambda, title) {
  ci <- calculate_ci(log_likelihoods, lambda_seq)
  plot(lambda_seq, log_likelihoods, type = "l", col = "blue", lwd = 2,
       main = title, xlab = expression(lambda), ylab = "Log-Verosimilitud")
  abline(h = ci$cutoff, col = "black", lty = 2)
  abline(v = ci$ci_lambda, col = "black", lty = 3)
  abline(v = optimal_lambda, col = "red", lty = 2)
}

# Aplicar y graficar para Masa Estelar
lambda_seq <- seq(-2, 2, 0.1)
log_likelihoods_masa <- result_masa$log_likelihoods
plot_log_likelihood(log_likelihoods_masa, lambda_seq, lambda_masa,
                    title = "Log-Verosimilitud para Masa Estelar")

# Aplicar y graficar para Luminosidad
log_likelihoods_luminosidad <- result_luminosidad$log_likelihoods
plot_log_likelihood(log_likelihoods_luminosidad, lambda_seq, lambda_luminosidad,
                    title = "Log-Verosimilitud para Luminosidad")

head(data)

## # A tibble: 6 × 14
##   bestobjid specobjid plate   mjd fiberid T_eff   Fe_H Log_g   snr flag 
##       <dbl>     <dbl> <dbl> <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl> <chr>
## 1   1.24e18   3.28e17   291 51928     321 5607. -0.595  4.22  40.2 nnnnn
## 2   1.24e18   3.28e17   291 51928     374 5551. -0.277  4.10  35.8 nnBnn
## 3   1.24e18   3.28e17   291 51928     426 5802. -0.785  4.20  34.5 nnnnn
## 4   1.24e18   3.28e17   291 51928     471 6206. -1.48   4.11  52.5 nnnnn
## 5   1.24e18   3.64e17   323 51615       2 4235. -1.23   4.30  30.9 nnngn
## 6   1.24e18   3.64e17   323 51615       5 4177. -1.50   4.48  30.4 nnngn
## # ℹ 4 more variables: Masa_estelar <dbl>, Luminosidad <dbl>,
## #   Masa_estelar_BoxCox <dbl>, Luminosidad_BoxCox <dbl>

Comparación de Distribuciones

Se comparan las distribuciones originales y transformadas de Masa y Luminosidad.

# Comparar distribuciones de Masa Estelar antes y después de Box-Cox
ggplot(data) +
  geom_density(aes(x = Masa_estelar), color = "blue", fill = "blue", alpha = 0.4) +
  geom_density(aes(x = Masa_estelar_BoxCox), color = "red", fill = "red", alpha = 0.4) +
  labs(
    title = "Distribución de Masa Estelar (Original vs Transformada)",
    x = "Masa Estelar",
    y = "Densidad"
  ) +
  scale_x_continuous(
    limits = c(min(c(data$Masa_estelar, data$Masa_estelar_BoxCox)),
               max(c(data$Masa_estelar, data$Masa_estelar_BoxCox))),
    breaks = scales::pretty_breaks(n = 10)
  ) +
  theme_minimal() +
  annotate("text", x = max(c(data$Masa_estelar, data$Masa_estelar_BoxCox)) - 0.5,
           y = max(density(data$Masa_estelar)$y) - 0.05, label = "Masa Estelar original", color = "blue", hjust = 1) +
  annotate("text", x = max(c(data$Masa_estelar, data$Masa_estelar_BoxCox)) - 0.5,
           y = max(density(data$Masa_estelar)$y) - 0.20, label = "Masa Estelar transformada", color = "red", hjust = 1)

# Comparar distribuciones de Luminosidad antes y después de Box-Cox
ggplot(data) +
  geom_density(aes(x = Luminosidad), color = "green", fill = "green", alpha = 0.4) +
  geom_density(aes(x = Luminosidad_BoxCox), color = "purple", fill = "purple", alpha = 0.4) +
  labs(
    title = "Distribución de Luminosidad (Original vs Transformada)",
    x = "Luminosidad",
    y = "Densidad"
  ) +
  scale_x_continuous(
    limits = c(-5, 15),
    breaks = scales::pretty_breaks(n = 10)
  ) +
  theme_minimal() +
  annotate("text", x = 14, y = 0.6, label = "Luminosidad original", color = "green", hjust = 1) +
  annotate("text", x = 14, y = 0.55, label = "Luminosidad transformada", color = "purple", hjust = 1)

## Warning: Removed 152 rows containing non-finite outside the scale range
## (`stat_density()`).

Comparación de Modelos Lineales

Se ajustan modelos de regresión lineal antes y después de la transformación para evaluar mejoras.

# Regresión antes de Box-Cox
modelo_original <- lm(Luminosidad ~ Masa_estelar, data = data)
summary(modelo_original)

## 
## Call:
## lm(formula = Luminosidad ~ Masa_estelar, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.209 -1.677 -0.118  1.301 47.790 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -11.9399     0.1756  -67.99   <2e-16 ***
## Masa_estelar  12.3589     0.1122  110.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.262 on 1998 degrees of freedom
## Multiple R-squared:  0.8587, Adjusted R-squared:  0.8586 
## F-statistic: 1.214e+04 on 1 and 1998 DF,  p-value: < 2.2e-16

# Regresión después de Box-Cox
modelo_transformado <- lm(Luminosidad_BoxCox ~ Masa_estelar_BoxCox, data = data)
summary(modelo_transformado)

## 
## Call:
## lm(formula = Luminosidad_BoxCox ~ Masa_estelar_BoxCox, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.183767 -0.048295  0.008508  0.055344  0.091020 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.091095   0.001444   -63.1   <2e-16 ***
## Masa_estelar_BoxCox  3.211409   0.003932   816.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06323 on 1998 degrees of freedom
## Multiple R-squared:  0.997,  Adjusted R-squared:  0.997 
## F-statistic: 6.669e+05 on 1 and 1998 DF,  p-value: < 2.2e-16

Resultados Finales

Se evalúan las predicciones realizadas a partir de los modelos ajustados.

plot_original <- ggplot(data, aes(x = Masa_estelar, y = Luminosidad)) +
  geom_point(alpha = 0.5, color = "blue") +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Regresión Lineal Antes de la Transformación",
    x = "Masa Estelar (Original)",
    y = "Luminosidad (Original)"
  ) +
  theme_minimal()

plot_transformado <- ggplot(data, aes(x = Masa_estelar_BoxCox, y = Luminosidad_BoxCox)) +
  geom_point(alpha = 0.5, color = "green") +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Regresión Lineal Después de la Transformación Box-Cox",
    x = "Masa Estelar (Box-Cox)",
    y = "Luminosidad (Box-Cox)"
  ) +
  theme_minimal()

summary_original <- summary(modelo_original)
cat("Modelo Antes de Box-Cox\n")

## Modelo Antes de Box-Cox

cat("R^2:", summary_original$r.squared, "\n")

## R^2: 0.8586526

cat("Error Estándar Residual (RMSE):", sqrt(mean(residuals(modelo_original)^2)), "\n\n")

## Error Estándar Residual (RMSE): 4.259427

summary_transformado <- summary(modelo_transformado)
cat("Modelo Después de Box-Cox\n")

## Modelo Después de Box-Cox

cat("R^2:", summary_transformado$r.squared, "\n")

## R^2: 0.997013

cat("Error Estándar Residual (RMSE):", sqrt(mean(residuals(modelo_transformado)^2)), "\n")

## Error Estándar Residual (RMSE): 0.06319512

print(plot_original)

## `geom_smooth()` using formula = 'y ~ x'

print(plot_transformado)

## `geom_smooth()` using formula = 'y ~ x'

## Validación con Nuevos Datos

En esta sección, se utilizan datos nuevos para evaluar el modelo desarrollado. El objetivo es predecir la luminosidad estelar en función de la masa estelar y validar las predicciones mediante el cálculo de residuos, \(R^2\), y MSE.

Cargar Nuevos Datos

Se cargan los nuevos datos desde un archivo .ods, los cuales contienen la temperatura efectiva (\(T_{\text{eff}}\)) para predicción.

file_path <- "/home/cesar/Documentos/Masas-a-predecir.ods" # Cambiar la ruta acorde a la ubicación en tu computadora
nuevos_datos <- read_ods(file_path, sheet = 1)
head(nuevos_datos)

## # A tibble: 6 × 10
##   bestobjid specobjid plate   mjd fiberid T_eff   Fe_H Log_g   snr flag 
##       <dbl>     <dbl> <dbl> <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl> <chr>
## 1   1.24e18   6.38e17   567 52252     282 4291. -0.319  4.02  36.8 nnnnn
## 2   1.24e18   9.50e17   844 52378     332 4428. -1.00   4.45  34.3 nnngn
## 3   1.24e18   9.50e17   844 52378     482 4450. -0.283  4.06  36.5 nnnnn
## 4   1.24e18   9.52e17   845 52381     613 4433. -0.474  4.31  34.4 nnnnn
## 5   1.24e18   9.55e17   848 52669      39 4444. -0.599  4.43  41.2 nnnnn
## 6   1.24e18   9.57e17   850 52338     318 4374. -1.15   4.42  48.2 nnngn

Aplicar Transformaciones Box-Cox

Se aplican las transformaciones Box-Cox usando los valores óptimos de \(\lambda\) obtenidos previamente.

inverse_boxcox <- function(y, lambda) {
  if (lambda == 0) {
    return(exp(y))
  } else {
    return((y * lambda + 1)^(1 / lambda))
  }
}

file_path <- "/home/cesar/Documentos/Masas-a-predecir.ods" 
nuevos_datos <- read_ods(file_path, sheet = 1)

if (!"T_eff" %in% colnames(nuevos_datos)) {
  stop("La columna 'T_eff' no está presente en los datos.")
}
nuevos_datos$T_eff <- as.numeric(nuevos_datos$T_eff)

nuevos_datos$Masa_estelar <- sapply(nuevos_datos$T_eff, calculate_mass)
nuevos_datos$Luminosidad <- mapply(calculate_luminosity, nuevos_datos$Masa_estelar, nuevos_datos$T_eff)

lambda_masa <- -0.7878788 
if (lambda_masa == 0) {
  nuevos_datos$Masa_estelar_BoxCox <- log(nuevos_datos$Masa_estelar)
} else {
  nuevos_datos$Masa_estelar_BoxCox <- (nuevos_datos$Masa_estelar^lambda_masa - 1) / lambda_masa
}

lambda_luminosidad <- -0.1818182
nuevos_datos$Luminosidad_Pred_Transformada <- predict(modelo_transformado, newdata = nuevos_datos)

nuevos_datos$Luminosidad_Pred <- inverse_boxcox(nuevos_datos$Luminosidad_Pred_Transformada, lambda_luminosidad)

head(nuevos_datos[c("T_eff", "Masa_estelar", "Masa_estelar_BoxCox", "Luminosidad", "Luminosidad_Pred")])

## # A tibble: 6 × 5
##   T_eff Masa_estelar Masa_estelar_BoxCox Luminosidad Luminosidad_Pred
##   <dbl>        <dbl>               <dbl>       <dbl>            <dbl>
## 1 4291.        0.551              -0.759       0.117            0.125
## 2 4428.        0.587              -0.661       0.147            0.156
## 3 4450.        0.593              -0.646       0.153            0.161
## 4 4433.        0.589              -0.658       0.148            0.157
## 5 4444.        0.592              -0.650       0.151            0.160
## 6 4374.        0.573              -0.699       0.135            0.143

Comparación de Predicciones

Se comparan las predicciones realizadas por el modelo con los valores calculados para validar el desempeño.

ggplot(nuevos_datos, aes(x = Luminosidad, y = Luminosidad_Pred)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  labs(
    title = "Comparación entre Luminosidad Calculada y Predicha",
    x = "Luminosidad Calculada",
    y = "Luminosidad Predicha"
  ) +
  theme_minimal()

### Métricas de Validación Se calculan \(R^2\) y el Error Cuadrático Medio (MSE) para evaluar el desempeño del modelo.

nuevos_datos$Luminosidad_Real <- nuevos_datos$Luminosidad
r_squared <- cor(nuevos_datos$Luminosidad_Real, nuevos_datos$Luminosidad_Pred)^2
mse <- mean((nuevos_datos$Luminosidad_Real - nuevos_datos$Luminosidad_Pred)^2)

cat("R²:", r_squared, "\n")

## R²: 0.9815023

cat("MSE:", mse, "\n")

## MSE: 4.294422

Análisis de Residuos

Se analiza el comportamiento de los residuos para identificar posibles sesgos en el modelo.

nuevos_datos$residuos <- nuevos_datos$Luminosidad_Real - nuevos_datos$Luminosidad_Pred

ggplot(nuevos_datos, aes(x = Luminosidad_Pred, y = residuos)) +
  geom_point(color = "darkorange", alpha = 0.7) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(
    title = "Gráfico de Residuos",
    x = "Luminosidad Predicha",
    y = "Residuos (Real - Predicha)"
  ) +
  theme_minimal()

Transformación Box-Cox para predecir luminosidad a partir de la masa estelar

Aguirre César, Bartolon Gustavo & Luján Uziel