library(data.table)
library(ggplot2)
library(caret)
library(jtools)
library(scales)
library(dplyr)
library(lattice)
library(MASS)
library(faraway)
library(cvTools)
library(lmtest)
library(faraway)
library(latex2exp)
library(cvTools)
library(Metrics)
rm(list = ls())

Pregunta 2

Cargado base de datos

df<- data.table(read.table('ew_f1_2.6598863068_ogle051153.02.txt', header = T, sep = " " , encoding = 'Latin-1'))

df<-df[,-3]

#Planteo inicial de datos#
    colnames(df)[1]<-"ti"
    colnames(df)[2]<-"yi"
    
    # Parametros Fijos
    f1= 2.6598863068
    # Pi ya se encuentra definido

Parte 2.a

Modelo Armónico

graph1 <- ggplot() +
  
  geom_point(df, mapping = aes(x=df$ti, y=df$yi, color='Valores reales')) +
  geom_line(df, mapping = aes(x=df$ti, y=df$yi), color='#FF36C2', method='loess',se=F, fullrange = T, alpha=0.3)+
  theme_minimal() +
  labs(x='ti',y='yi',color=NULL) +
  theme(legend.position = 'bottom')
## Warning in geom_line(df, mapping = aes(x = df$ti, y = df$yi), color =
## "#FF36C2", : Ignoring unknown parameters: `method`, `se`, and `fullrange`
graph1

Este gráfico presenta el ajuste de la línea de luz original. Incorpora los valores reales y los ajustes pertinentes para denotar el comportamiento de los datos.

Modelo Armónico Tipo 1

Implementación propia

# Define the value of f1 and k

k <- 6

# Create a data frame to store the harmonic terms
harmonics <- data.frame(matrix(0, nrow = nrow(df), ncol = 2 * k))

# Generate the harmonic terms
for (j in 1:k) {
  harmonics[, (2 * j - 1)] <- sin(2 * pi * f1 * j * df$ti)
  harmonics[, (2 * j)] <- cos(2 * pi * f1 * j * df$ti)
}

# Combine the harmonic terms with the other variables
df_harmonics <- cbind(df, harmonics)


# Perform nonlinear regression
reg1 <- nls(yi ~ b0 + b1 * ti + rowSums(harmonics * rep(c(a, b), each = nrow(df))),
            data = df, start = c(a = 1, b = 1, b0 = 0, b1 = 0))

sum01<-summary(reg1)
sum01
## 
## Formula: yi ~ b0 + b1 * ti + rowSums(harmonics * rep(c(a, b), each = nrow(df)))
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## a  1.533e-02  2.401e-03   6.388 5.94e-10 ***
## b  1.522e-02  2.562e-03   5.944 7.31e-09 ***
## b0 4.973e+00  2.404e+01   0.207    0.836    
## b1 4.963e-06  9.807e-06   0.506    0.613    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07464 on 319 degrees of freedom
## 
## Number of iterations to convergence: 1 
## Achieved convergence tolerance: 1.933e-07

Este modelo no se considera más adelante, puesto que considera un funcionamiento en el cual solo se incorpora y analiza el valor de un K a la vez. Por ende no cumple con el requerimiento de mostrar la totalidad de ciclos respecto al valor máximo del K. No obstante, permite analizar correctamente el ajuste de cada una de las iteraciones incluidas.

Grafico modelo Armónico Tipo 1

graph1 <- ggplot() +
  
    geom_point(df, mapping = aes(x=df$ti, y=df$yi, color='Valores Reales')) +
  geom_line(df, mapping = aes(x=df$ti, y=predict(reg1)), color='#36FFA4', method='loess',se=F, fullrange = T) +
  geom_line(df, mapping = aes(x=df$ti, y=df$yi), color='#FF36C2', method='loess',se=F, fullrange = T, alpha=0.3)+
  
  theme_minimal() +
  
  labs(x='ti',y='yi',color=NULL) +
  
  theme(legend.position = 'bottom') +
  
  scale_y_continuous(labels = number_format(scale = 1))

graph1

Calculo MSE Modelo 1

mean((df$yi-predict(reg1))^2)
## [1] 0.005501552

La idea de aplicación de este modelo se basa en la idea de no linealidad. Por ende, incorpora funciones y analiza su comportamiento dependiendo del valor K asignado arbitrariamente, pero carece de la incorporación o asignación de los demás Y asociados.

Modelo Armónico Tipo 2

Aplicación y modificación de Ayudantía

## Determinación de K y generación de variables adicionales
k<-1

for (j in 1:k) {
  componente_a <-  sin(2 * pi * f1 *j* df$ti) 
  componente_b <- cos(2 * pi * f1 *j* df$ti)
  col_name_a <- paste0("componente_a")
  col_name_b <- paste0("componente_b")
  df$col_name_a <- componente_a
  df$col_name_b <- componente_b
}
# Implementación de la regresión
f2=formula(yi~ ti + col_name_a + col_name_b)


reg2 = lm(formula=f2, data=df)

pred <- predict(reg2)

sum02<-summary((reg2))
sum02
## 
## Call:
## lm(formula = f2, data = df)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.065736 -0.017436 -0.000261  0.018424  0.074491 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.841e+01  9.013e+00   2.043   0.0419 *  
## ti          -5.223e-07  3.677e-06  -0.142   0.8871    
## col_name_a   6.500e-02  2.173e-03  29.916   <2e-16 ***
## col_name_b   8.884e-02  2.242e-03  39.629   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02802 on 319 degrees of freedom
## Multiple R-squared:  0.8856, Adjusted R-squared:  0.8845 
## F-statistic: 822.8 on 3 and 319 DF,  p-value: < 2.2e-16

Grafico Modelo tipo 2

graph1 <- ggplot() +
  
    geom_point(df, mapping = aes(x=df$ti, y=df$yi, color='Valores Reales')) +
  geom_line(df, mapping = aes(x=df$ti, y=predict(reg2)), color='#36FFA4') +
  geom_line(df, mapping = aes(x=df$ti, y=df$yi), color='#FF36C2', alpha=0.3)+
  
  theme_minimal() +
  
  labs(x='ti',y='yi',color=NULL) +
  
  theme(legend.position = 'bottom') +
  
  scale_y_continuous(labels = number_format(scale = 1))

graph1

Calculo MSE Modelo 2

mean((df$yi-predict(reg2))^2)
## [1] 0.0007751282

Este modelo, a diferencia del modelo anterior si considera la interacción de las variables dependiendo del k. Por ende realiza una asignación conjunta para realizar el análisis dependiendo del número de valores asignados. No obstante, parece que tampoco cumple con los valores y requerimientos solicitados por la pregunta.

Modelo Armónico Tipo

Aplicación y modificación de Ayudantía

# Define the number of terms and f1 value
k <- 4

# Initialize an empty formula string
formula <- "df$yi ~ df$ti"

# Iterate over the terms and add them to the formula string
for (i in 1:k) {
  # Generate the sine and cosine terms for each i
  sin_term <- paste("sin(2 * pi * f1 *", i, "* df$ti)", sep = "")
  cos_term <- paste("cos(2 * pi * f1 *", i, "* df$ti)", sep = "")
  
  # Add the sine and cosine terms to the formula string
  formula <- paste(formula, "+", sin_term, "+", cos_term, sep = " ")
  
  
  df[, paste0("sin_term", i)] <- eval(parse(text = sin_term))
  df[, paste0("cos_term", i)] <- eval(parse(text = cos_term))
  
}

# Create the final formula
f3 <- formula(yi~., data=df)
reg3 = lm(formula=f3, data=df)

pred <- predict(reg3)

sum03<-summary((reg3))
sum03
## 
## Call:
## lm(formula = f3, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.04670 -0.01065 -0.00139  0.01333  0.06640 
## 
## Coefficients: (2 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.790e+01  5.863e+00   3.052  0.00246 ** 
## ti          -3.111e-07  2.392e-06  -0.130  0.89662    
## col_name_a   6.418e-02  1.416e-03  45.336  < 2e-16 ***
## col_name_b   8.995e-02  1.465e-03  61.397  < 2e-16 ***
## sin_term1           NA         NA      NA       NA    
## cos_term1           NA         NA      NA       NA    
## sin_term2    2.594e-02  1.404e-03  18.474  < 2e-16 ***
## cos_term2    7.649e-03  1.476e-03   5.181 3.97e-07 ***
## sin_term3    1.217e-02  1.455e-03   8.359 2.10e-15 ***
## cos_term3   -4.414e-03  1.428e-03  -3.091  0.00218 ** 
## sin_term4    3.321e-03  1.442e-03   2.303  0.02194 *  
## cos_term4   -2.058e-03  1.440e-03  -1.429  0.15399    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01817 on 313 degrees of freedom
## Multiple R-squared:  0.9527, Adjusted R-squared:  0.9514 
## F-statistic: 701.1 on 9 and 313 DF,  p-value: < 2.2e-16
graph1 <- ggplot() +
  
    geom_point(df, mapping = aes(x=df$ti, y=df$yi, color='Valores Reales')) +
  geom_line(df, mapping = aes(x=df$ti, y=predict(reg3)), color='#36FFA4') +
  geom_line(df, mapping = aes(x=df$ti, y=df$yi), color='#FF36C2', alpha=0.3)+
  
  theme_minimal() +
  
  labs(x='ti',y='yi',color=NULL) +
  
  theme(legend.position = 'bottom') +
  
  scale_y_continuous(labels = number_format(scale = 1))

graph1

Calculo MSE Modelo 3

mean((df$yi-predict(reg3))^2)
## [1] 0.0003200887

Este modelo final, incorpora la esencia del planteo armónico. Se enfoca en considerar las distintas instancias de la no linealidad e incorpora los parámetros lineales estimados para cada uno dentro de la estimación. Por lo cual, a diferencia de los modelos anteriores, minimiza el error cuadrático medio con un K=4, a diferencia de los primeros dos modelos que lo encontraban en un K=1.

Parte 2.b

Como se ha presentado en la parte anterior, se puede apreciar que en los gráficos el ajuste del modelo 3 supera a los dos primeros en la línea ajustada para cada parámetro. No obstante, cabe mencionar que todos los modelos en esencia cuentan con un sobre ajuste de los datos.

Por otro lado, la razón principal de la elección del K=4, junto al modelo 3 es que este modelo minimiza el error cuadrático medio, junto a lo que el K va a determinar el valor mínimo que puede tomar sin entrar en un sobre ajuste excesivo de los datos.

Adicionalmente, se considera que al aumentar el K valores sobre 4, esto produce que el sobre ajuste quite peso a variables importantes para el análisis y afecta en el modelo en general, ya que las variables incorporadas carecen de significancia para el modelo y se intenta minimizar el nivel de variables omitidas.

Parte 2.c

Implementación K-folds CV

Modelo 2 y Modelo 3

data.table(Regresión=c('Modelo 03','Modelo 02'),
           
           `R2`=c(sum03$r.squared, sum02$r.squared),
           
           `R2 ajustado`=c(sum03$adj.r.squared, sum02$adj.r.squared),
           
           `Mean Absolute Error` = c(MAE(reg3$fitted.values, df$yi), 
                                     MAE(reg2$fitted.values, df$yi)),
           
           `Root Mean Squared Error`=c(RMSE(reg3$fitted.values, df$yi),
                                      RMSE(reg2$fitted.values, df$yi)))

Implementación K-Folds Modelo 3

set.seed(12345) ## setear una semilla

training.samples <- createDataPartition(df$yi,p = 0.8, list = FALSE)

train.data  <- df[training.samples, ]

test.data <- df[-training.samples, ]

model <- lm(formula = f3, data =  train.data)

predictions01 <- predict(model,train.data)
predictions02 <- predict(model,test.data)

data.table(' ' = c('Dentro de muestra','Fuera de muestra'),

           RMSE = c(RMSE(predictions01, train.data$yi),
                    RMSE(predictions02, test.data$yi)),
           MAE = c(MAE(predictions01, train.data$yi),
                   MAE(predictions02, test.data$yi)))

Por lo tanto, parece que el modelo es posee un bajo nivel de explicación associado a las variables fuera de la muestra en el proceso de entrenamiento.

Implementación LOO CV

LOO CV Modelo 3

set.seed(123)## setear una semilla
setupLOO <- trainControl(method = "LOOCV")

predLOO <- train(f3,data=df[1:323],method="lm",trControl= setupLOO)
print(predLOO)
## Linear Regression 
## 
## 323 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation 
## Summary of sample sizes: 322, 322, 322, 322, 322, 322, ... 
## Resampling results:
## 
##   RMSE        Rsquared   MAE       
##   0.01845812  0.9497015  0.01466882
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Como presentamos en este análisis, el modelo se somete a dos formas de validación cruzada para la determinación de su error de predicción. En este caso va a corresponder a un 0.01845812 con el método de LOO y a un 0.02167478 con el método de K-Folds considerando los valores fuera de la muestra.

Por estas razones el modelo muestra un buen desempeño en las pruebas actuales y es un candidato idóneo para el desarrollo de consideraciones dentro de la muestra. No obstante, es necesario someterlo a nuevas validaciones y formas para poder justificar la elección correcta del K y la formulación escogida.

Pregunta 3

Cargado base de datos

df<- data.table(read.table('ew_f1_2.6598863068_ogle051153.02.txt', header = T, sep = " " , encoding = 'Latin-1'))

df<-df[,-3]

#Planteo inicial de datos#
    colnames(df)[1]<-"ti"
    colnames(df)[2]<-"yi"
    

    # Pi ya se encuentra definido

Sección 3.a

Transformación solicitada implementada

Para la consideración del tiempo t0 se va a considerar el mínimo correspondiente a la variable del tiempo. Esto garantiza que sea el tiempo inicial en un instante propio del tiempo. No obstante, este resultado se va a modificar en el caso de utilizar directamente el valor de ti=2450458, cambiando la concavidad del problema. Adicionalmente, si se toma en consideración el tiempo inicial con la sustracción de una unidad para asegurar un instante de tiempo fuera de la base, este modelo también se modifica. Se incorporan comentarios con los códigos para realizar las pruebas pertinentes.

  # Definir los parámetros t0 y f1
  t0 <- min(df$ti)  # Valor de referencia para t0

              #t0 <- 2450458  # Valor directo
            
              #t0 <- min(df$ti)-1  # Valor fuera de la base

    # Parametros Fijos
    f1= 2.6598863068

# Calcular zi
df$zi <- (df$ti - t0) %% (1 / f1)

# Graficar zi vs yi

graph1 <- ggplot() +
  
    geom_point(df, colour = "#36DAFF", size = 1, mapping = aes(x=df$zi, y=df$yi)) +
  
  theme_minimal() +
  
  labs(x='zi',y='yi', title =  "Gráfico de zi v/s yi") +
  
  theme(legend.position = 'bottom') +
  
  scale_y_continuous(labels = number_format(scale = 1))

graph1

Sección 3.b

Calculo de ajuste polinomial

# Ajuste de polinomio a los datos (zi, yi)
grados <- 1:9  # Rango de grados del polinomio a probar
ajustes <- list()  # Lista para almacenar los ajustes

# Iterar sobre los diferentes grados del polinomio
for (grado in grados) {
  polinomio <- poly(df$zi, degree = grado, raw = TRUE)  # Generar polinomio de grado 'grado'
  ajuste <- lm(df$yi ~ polinomio)  # Ajustar modelo lineal utilizando el polinomio
  
  # Almacenar el ajuste en la lista
  ajustes[[as.character(grado)]] <- ajuste
}

# Calcular los errores cuadrados medios (MSE) de los ajustes
mse <- sapply(ajustes, function(x) mean(x$residuals^2))

# Identificar el menor grado del polinomio que ajusta a los datos
menor_grado <- which.min(mse)

menor_grado
## 9 
## 9

Implementación de modelo polinomial

# Generar el polinomio de grado 10

pol_10 <- poly(df$zi, degree = menor_grado, raw = TRUE)

f4<-yi ~ pol_10


# Ajustar el modelo de regresión polinomial con grado 10
reg_pol <- lm(f4, data=df)

df$pred<-predict(reg_pol)

# Imprimir el resumen del modelo
sum04<-summary(reg_pol)
sum04
## 
## Call:
## lm(formula = f4, data = df)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.046889 -0.011062 -0.001053  0.013358  0.066347 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  1.726e+01  8.093e-03 2132.913   <2e-16 ***
## pol_101      3.354e+00  1.321e+00    2.539   0.0116 *  
## pol_102     -1.019e+02  6.837e+01   -1.491   0.1369    
## pol_103      1.680e+02  1.595e+03    0.105   0.9162    
## pol_104      1.375e+04  2.001e+04    0.687   0.4925    
## pol_105     -1.690e+05  1.465e+05   -1.154   0.2493    
## pol_106      9.337e+05  6.439e+05    1.450   0.1480    
## pol_107     -2.755e+06  1.674e+06   -1.646   0.1007    
## pol_108      4.213e+06  2.368e+06    1.779   0.0762 .  
## pol_109     -2.626e+06  1.405e+06   -1.869   0.0625 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01817 on 313 degrees of freedom
## Multiple R-squared:  0.9528, Adjusted R-squared:  0.9514 
## F-statistic: 701.6 on 9 and 313 DF,  p-value: < 2.2e-16

Grafico de regresión polinomial

ggplot() +
  geom_segment(df, mapping = aes(x=df$zi, xend=df$zi, y=df$yi, yend=df$pred)) +
  geom_point(df, mapping = aes(x=df$zi, y=df$yi, color='Valores reales')) +
  geom_line(df, mapping = aes(x=df$zi, y=df$pred, color='Valores predichos'), method='loess',se=T, fullrange = T, size=1) +
  theme_minimal() +
  labs(x='Position',y='Height',color=NULL, title = "Modelo 2") +
  theme(legend.position = 'bottom') +
  scale_y_continuous(labels = number_format(scale = 1))

El proceso de selección del K se realizó gracias a la consideración del error cuadrático medio. En este modelo se escoge e implementa un K=9, no obstante, se ha determinado que para minimizar el sobre ajuste de los datos y tener un modelo más generalizable es necesario tomar un K entre el valor de 2 y el valor de 4.

No obstante, para buscar generar la comparación más cercana con el modelo anterior, se ha escogido el valor de 9, que ajusta los valores sin llegar a un sobre ajuste excesivo, pasando por el centro de toda la muestra, minimizando el error cuadrático bajo las consideraciones antes expuestas.

Sección 3.c

K-Folds Modelo Anterior

set.seed(12345) ## setear una semilla

training.samples <- createDataPartition(df$yi,p = 0.8, list = FALSE)

train.data  <- df[training.samples, ]

test.data <- df[-training.samples, ]

model <- lm(f3, data = train.data)

predictions01 <- predict(model,train.data)
predictions02 <- predict(model,test.data)

data.table(' ' = c('Dentro de muestra','Fuera de muestra'),

           RMSE = c(RMSE(predictions01, train.data$yi),
                    RMSE(predictions02, test.data$yi)),
           MAE = c(MAE(predictions01, train.data$yi),
                   MAE(predictions02, test.data$yi)))

K-Folds Modelo Adaptado

set.seed(12345) ## setear una semilla

training.samples <- createDataPartition(df$yi,p = 0.8, list = FALSE)

train.data  <- df[training.samples, ]

test.data <- df[-training.samples, ]

model <- lm(yi~ zi + zi^2 + zi^3 + zi^4 + zi^5 + zi^6, data = train.data)

predictions01 <- predict(model,train.data)
predictions02 <- predict(model,test.data)

data.table(' ' = c('Dentro de muestra','Fuera de muestra'),

           RMSE = c(RMSE(predictions01, train.data$yi),
                    RMSE(predictions02, test.data$yi)),
           MAE = c(MAE(predictions01, train.data$yi),
                   MAE(predictions02, test.data$yi)))

Consideración de CV conjunta

data.table(Regresión=c('Modelo 03','Modelo 04'),
           
           `R2`=c(sum03$r.squared, sum04$r.squared),
           
           `R2 ajustado`=c(sum03$adj.r.squared, sum04$adj.r.squared),
           
           `Mean Absolute Error` = c(MAE(reg3$fitted.values, df$yi), 
                                     MAE(reg_pol$fitted.values, df$yi)),
           
           `Root Mean Squared Error`=c(RMSE(reg3$fitted.values, df$yi),
                                      RMSE(reg_pol$fitted.values, df$yi)))

En esta sección se presentan todos los métodos de evaluación implementados para la determinación del modelo. Dadas las características del modelo 3 y 4, es esperable que exista un mejor ajuste con el modelo 3, lo cual es el resultado obtenido del análisis. No obstante, la diferencia entre el modelo transformado y el modelo básico es mínima. Por lo cual, es mejor presentar la diferencia con un modelo menos ajustado en la parte 3. Ya que esto va a dar mejor punto de comparación para la determinación final.

Pregunta 4

Sección 4.a

Cargado base de datos

rm(list = ls())
df<- data.table(read.table('sys.txt', header = T, sep = " " , encoding = 'Latin-1'))

    colnames(df)[1]<-"sys1"

Modelo de Componentes Principales

library(pls)
## 
## Attaching package: 'pls'
## The following object is masked from 'package:caret':
## 
##     R2
## The following object is masked from 'package:stats':
## 
##     loadings
library(psych)
## 
## Attaching package: 'psych'
## The following object is masked from 'package:faraway':
## 
##     logit
## The following objects are masked from 'package:scales':
## 
##     alpha, rescale
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
#Split the data

  idx<-floor(0.8*nrow(df))

#Set the partition

  train_ind<-sample(seq_len(nrow(df)), size=idx)
    train<-df[train_ind,]
    test<-df[-train_ind,]
    pc<-df[,-1]

Modelo

pc.fit<-prcomp(~., data=train)

screeplot(pc.fit, type="l", main="Sreeplot for SYS features" )

# Estimar el modelo de regresión en componentes principales
modelo_pls <- plsr(sys1 ~ .,ncomp=1, data=df, validation="CV")  # ncomp es el número de componentes principales a considerar

# Resumen del modelo
sum00<-summary(modelo_pls)
## Data:    X dimension: 1643 2658 
##  Y dimension: 1643 1
## Fit method: kernelpls
## Number of components considered: 1
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps
## CV           20.82    21.16
## adjCV        20.82    21.12
## 
## TRAINING: % variance explained
##       1 comps
## X       2.352
## sys1    4.464
sum00
## NULL
# Predicción
pls_predict<-predict(modelo_pls)


# Calculo de Error cuadratico medio
rmse_value <- rmse(actual=df$sys1, predicted=pls_predict)


# R-squared (R2) calculation
y_mean <- mean(df$sys1)
ss_total <- sum((df$sys1 - y_mean)^2)
ss_residual <- sum((df$sys1 - pls_predict)^2)
r_squared <- 1 - (ss_residual / ss_total)

# Adjusted R-squared (R2 adjusted) calculation
n <- nrow(train)
p <- ncol(train) - 1  # Number of predictors (excluding the intercept)

r_squared_adjusted <- (1-(ss_residual / ss_total)*((n-p)/(n-1)))

# Print RMSE, R-squared, and adjusted R-squared
print(paste("RMSE:", rmse_value))
## [1] "RMSE: 20.3387868916449"
print(paste("R-squared:", r_squared))
## [1] "R-squared: 0.0446448036865995"
print(paste("Adjusted R-squared:", r_squared_adjusted))
## [1] "Adjusted R-squared: 1.97791118343123"

En este modelo se aplican las conspiraciones indagatorias antes realizadas para poder determinar el número de componentes. En esta formulación se implementa un modelo de componentes que minimiza el error y toma en cuenta los valores de prueba para el modelo final.

Sección 4.b

Implementación de modelo Ridge

library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 4.1-7
#define response variable
y <- data.matrix(df[,1])

#define matrix of predictor variables
x <- data.matrix(df[,-1])
#fit ridge regression model
model <- glmnet(x, y, alpha = 0)

#view summary of model
summary(model)
##           Length Class     Mode   
## a0           100 -none-    numeric
## beta      265800 dgCMatrix S4     
## df           100 -none-    numeric
## dim            2 -none-    numeric
## lambda       100 -none-    numeric
## dev.ratio    100 -none-    numeric
## nulldev        1 -none-    numeric
## npasses        1 -none-    numeric
## jerr           1 -none-    numeric
## offset         1 -none-    logical
## call           4 -none-    call   
## nobs           1 -none-    numeric

The optimal lamda

#perform k-fold cross-validation to find optimal lambda value
cv_model <- cv.glmnet(x, y, alpha = 0)

#find optimal lambda value that minimizes test MSE
best_lambda <- cv_model$lambda.min
best_lambda
## [1] 1484.406

Plot considerando los lambas

#produce plot of test MSE by lambda value
plot(cv_model)

Fit con el mejor lambda

#find coefficients of best model
best_model <- glmnet(x, y, alpha = 0, lambda = best_lambda)
sum01<-summary(best_model)
sum01
##           Length Class     Mode   
## a0           1   -none-    numeric
## beta      2658   dgCMatrix S4     
## df           1   -none-    numeric
## dim          2   -none-    numeric
## lambda       1   -none-    numeric
## dev.ratio    1   -none-    numeric
## nulldev      1   -none-    numeric
## npasses      1   -none-    numeric
## jerr         1   -none-    numeric
## offset       1   -none-    logical
## call         5   -none-    call   
## nobs         1   -none-    numeric

#produce Ridge trace plot
plot(model, xvar = "lambda")
abline(h=0,
       col = "yellow",            # Modify color
       lty = "dashed",         # Modify line type
       lwd = 3)

**

#use fitted best model to make predictions
y_predicted <- predict(model, s = best_lambda, newx = x)

#find SST and SSE
sst <- sum((y - mean(y))^2)
sse <- sum((y-y_predicted)^2)

#find R-Squared
R2_ridge <- 1 - sse/sst
R2_ridge
## [1] 0

El modelo de Ridge presenta todos los cálculos asociados a su implementación y visualización. Por los resultados obtenidos se presenta una conclusión improbable. Ya que en esta formulación se muestra como la mayoría de las variables tienden a cero, esto produce que la conclusión final sea que la presión sistólica no pareciera tener tanta determinación a partir de una componente genética. Por lo cual, la alternativa es que el modelo puede mostrar ciertas variables influyentes, pero existe otra información que no está observada por la determinación de esta variable.

Sección 4.c

# R-squared (R2) calculation
y_mean1 <- mean(y)
ss_total1 <- sum((y - y_mean1)^2)
ss_residual1 <- sum((y - y_predicted)^2)
r_squared1 <- 1 - (ss_residual1 / ss_total1)

# Adjusted R-squared (R2 adjusted) calculation
n1 <- nrow(df)
p1 <- ncol(df) - 1  # Number of predictors (excluding the intercept)

adj_R2_ridge <- (1-(ss_residual / ss_total)*((n1-p1)/(n1-1)))

rmse_ridge <- rmse(actual=y, predicted=y_predicted)


# Print RMSE, R-squared, and adjusted R-squared
print(paste("RMSE:", rmse_ridge))
## [1] "RMSE: 20.8085879664976"
print(paste("R-squared:", r_squared1))
## [1] "R-squared: 0"
print(paste("Adjusted R-squared:", adj_R2_ridge))
## [1] "Adjusted R-squared: 1.59055147640566"
data.table(Regresión=c('Modelo CP','Modelo Ridge'),
           
           `R2`=c(r_squared, r_squared1),
           
           `R2 ajustado`=c(r_squared_adjusted, adj_R2_ridge),
           
           `Root Mean Squared Error`=c(rmse_value, rmse_ridge))

Dada la comparación cruzada extraída de las dos implementaciones, podemos concluir que el modelo de componentes principales muestra un mejor ajuste que el modelo de Ridge, en la consideración de poder explicar mejor los datos. Por otro lado, el modelo de Ridge nos permite apreciar la importancia de algunas variables en el análisis final, puesto que el modelo en sí penaliza ciertos valores y los acerca a cero.