library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ ggplot2 3.4.3 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
ENSANUT_2006 <- read_excel("C:/Users/Asus ZenBook/Downloads/ENSANUT 2006.xlsx",
sheet = "ENSANUT 2006")
ENSANUT_2006 <- mutate_all(ENSANUT_2006, ~ifelse(is.na(.), "No contestó", .))
Regresión lineal
# Regresión lineal
modelo <- lm(bmi ~ liquidos + miscelaneos + verduras + carnes + lacteos + cereales + azucar, data = ENSANUT_2006)
summary(modelo)
##
## Call:
## lm(formula = bmi ~ liquidos + miscelaneos + verduras + carnes +
## lacteos + cereales + azucar, data = ENSANUT_2006)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.742 -0.174 2.359 4.473 199.851
##
## Coefficients: (6 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.82948 0.20285 58.316 < 2e-16 ***
## liquidos2 0.18928 0.24528 0.772 0.44030
## liquidosNo contestó 1.91255 0.21246 9.002 < 2e-16 ***
## miscelaneos2 0.28029 0.24192 1.159 0.24664
## miscelaneosNo contestó NA NA NA NA
## verduras2 -0.08514 0.25942 -0.328 0.74276
## verdurasNo contestó NA NA NA NA
## carnes2 -0.30089 0.24222 -1.242 0.21418
## carnesNo contestó NA NA NA NA
## lacteos2 0.68866 0.24155 2.851 0.00436 **
## lacteosNo contestó NA NA NA NA
## cereales2 -0.59828 0.29405 -2.035 0.04190 *
## cerealesNo contestó NA NA NA NA
## azucar2 -0.57725 0.22174 -2.603 0.00924 **
## azucar9 -0.23549 0.80470 -0.293 0.76980
## azucarNo contestó NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.398 on 23995 degrees of freedom
## Multiple R-squared: 0.01215, Adjusted R-squared: 0.01178
## F-statistic: 32.78 on 9 and 23995 DF, p-value: < 2.2e-16
Regresión multiple
# Regresión multiple
multi_model <- lm(bmi ~ liquidos + cereales + verduras + carnes + lacteos + miscelaneos + azucar, data = ENSANUT_2006)
summary(multi_model)
##
## Call:
## lm(formula = bmi ~ liquidos + cereales + verduras + carnes +
## lacteos + miscelaneos + azucar, data = ENSANUT_2006)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.742 -0.174 2.359 4.473 199.851
##
## Coefficients: (6 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.82948 0.20285 58.316 < 2e-16 ***
## liquidos2 0.18928 0.24528 0.772 0.44030
## liquidosNo contestó 1.91255 0.21246 9.002 < 2e-16 ***
## cereales2 -0.59828 0.29405 -2.035 0.04190 *
## cerealesNo contestó NA NA NA NA
## verduras2 -0.08514 0.25942 -0.328 0.74276
## verdurasNo contestó NA NA NA NA
## carnes2 -0.30089 0.24222 -1.242 0.21418
## carnesNo contestó NA NA NA NA
## lacteos2 0.68866 0.24155 2.851 0.00436 **
## lacteosNo contestó NA NA NA NA
## miscelaneos2 0.28029 0.24192 1.159 0.24664
## miscelaneosNo contestó NA NA NA NA
## azucar2 -0.57725 0.22174 -2.603 0.00924 **
## azucar9 -0.23549 0.80470 -0.293 0.76980
## azucarNo contestó NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.398 on 23995 degrees of freedom
## Multiple R-squared: 0.01215, Adjusted R-squared: 0.01178
## F-statistic: 32.78 on 9 and 23995 DF, p-value: < 2.2e-16
summary_result_multi <- summary(multi_model)
p_values_multi <- summary_result_multi$coefficients[, "Pr(>|t|)"]
alpha_multi <- 0.05
for (i in seq_along(p_values_multi)) {
if (p_values_multi[i] < alpha_multi) {
cat(names(p_values_multi)[i], "is statistically significant (p-value <", alpha_multi, ")\n")
} else {
cat(names(p_values_multi)[i], "is not statistically significant (p-value >=", alpha_multi, ")\n")
}
}
## (Intercept) is statistically significant (p-value < 0.05 )
## liquidos2 is not statistically significant (p-value >= 0.05 )
## liquidosNo contestó is statistically significant (p-value < 0.05 )
## cereales2 is statistically significant (p-value < 0.05 )
## verduras2 is not statistically significant (p-value >= 0.05 )
## carnes2 is not statistically significant (p-value >= 0.05 )
## lacteos2 is statistically significant (p-value < 0.05 )
## miscelaneos2 is not statistically significant (p-value >= 0.05 )
## azucar2 is statistically significant (p-value < 0.05 )
## azucar9 is not statistically significant (p-value >= 0.05 )
Matriz de correlación
# Matriz de correlación
cor_matrix <- cor(ENSANUT_2006[, c("edad_anos", "bmi", "sexo")])
print(cor_matrix)
## edad_anos bmi sexo
## edad_anos 1.000000000 0.15261876 0.008228689
## bmi 0.152618759 1.00000000 -0.010481616
## sexo 0.008228689 -0.01048162 1.000000000
interpret_correlation <- function(cor_value) {
if (cor_value > 0.8) {
return("Muy Fuerte")
} else if (cor_value > 0.6) {
return("Fuerte")
} else if (cor_value > 0.4) {
return("Moderado")
} else if (cor_value > 0.2) {
return("Debil")
} else {
return("Muy debil")
}
}
interpreted_matrix <- apply(cor_matrix, c(1, 2), interpret_correlation)
print(interpreted_matrix)
## edad_anos bmi sexo
## edad_anos "Muy Fuerte" "Muy debil" "Muy debil"
## bmi "Muy debil" "Muy Fuerte" "Muy debil"
## sexo "Muy debil" "Muy debil" "Muy Fuerte"
heatmap(cor_matrix,
annot = interpreted_matrix, # show the interpretations in each cell
cmap = "coolwarm", # color map
main = "Correlation Matrix with Interpretations")
## Warning in plot.window(...): "annot" is not a graphical parameter
## Warning in plot.window(...): "cmap" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "annot" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "cmap" is not a graphical parameter
## Warning in title(...): "annot" is not a graphical parameter
## Warning in title(...): "cmap" is not a graphical parameter
Anova
modelo <- lm(bmi ~ liquidos + cereales + verduras + carnes + miscelaneos + lacteos + azucar + edad_anos, data = ENSANUT_2006)
anova_resultados <- anova(modelo)
print(anova_resultados)
## Analysis of Variance Table
##
## Response: bmi
## Df Sum Sq Mean Sq F value Pr(>F)
## liquidos 2 19281 9640.3 138.4141 < 2.2e-16 ***
## cereales 1 271 271.1 3.8918 0.048534 *
## verduras 1 1 1.3 0.0190 0.890291
## carnes 1 45 45.3 0.6501 0.420078
## miscelaneos 1 160 159.8 2.2945 0.129847
## lacteos 1 567 567.3 8.1450 0.004322 **
## azucar 2 482 241.2 3.4637 0.031329 *
## edad_anos 1 21256 21256.3 305.1954 < 2.2e-16 ***
## Residuals 23994 1671138 69.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Modelo Predictivo
print(paste("Modelo Predictivo - Alimentos"))
## [1] "Modelo Predictivo - Alimentos"
df_cleaned <- ENSANUT_2006[, c("bmi", "edad_anos", "liquidos", "cereales", "verduras", "carnes", "miscelaneos", "lacteos", "azucar")]
df_cleaned <- df_cleaned[complete.cases(df_cleaned), ]
set.seed(123)
split_index <- sample(1:nrow(df_cleaned), 0.8 * nrow(df_cleaned))
train_data <- df_cleaned[split_index, ]
test_data <- df_cleaned[-split_index, ]
lm_model <- lm(bmi ~ liquidos + cereales + verduras + carnes + miscelaneos + lacteos + azucar + edad_anos, data = train_data)
predictions <- predict(lm_model, newdata = test_data)
mse <- mean((test_data$bmi - predictions)^2)
r_squared <- summary(lm_model)$r.squared
adjusted_r_squared <- summary(lm_model)$adj.r.squared
mae <- mean(abs(test_data$bmi - predictions))
rmse <- sqrt(mse)
print(paste("Model Evaluation Metrics:"))
## [1] "Model Evaluation Metrics:"
print(paste("Mean Squared Error (MSE): ", mse))
## [1] "Mean Squared Error (MSE): 71.5174475137006"
print(paste("Mean Absolute Error (MAE): ", mae))
## [1] "Mean Absolute Error (MAE): 5.9651698842639"
print(paste("Root Mean Squared Error (RMSE): ", rmse))
## [1] "Root Mean Squared Error (RMSE): 8.45679889282586"
print(paste("R-squared: ", r_squared))
## [1] "R-squared: 0.0237437950273371"
print(paste("Adjusted R-squared: ", adjusted_r_squared))
## [1] "Adjusted R-squared: 0.023235142807792"
if (mse == 0) {
print("The model perfectly predicts the 'bmi' variable.")
} else {
print(paste("En promedio, las predicciones del modelo están equivocadas por aproximadamente", round(sqrt(mse), 2), "unidades en la escala de 'bmi'"))
}
## [1] "En promedio, las predicciones del modelo están equivocadas por aproximadamente 8.46 unidades en la escala de 'bmi'"
if (rmse < 10) {
cat("El RMSE es bajo, sugiriendo que el modelo tiene un buen ajuste.\n")
} else {
cat("El RMSE es alto, indicando que el modelo puede tener limitaciones.\n")
}
## El RMSE es bajo, sugiriendo que el modelo tiene un buen ajuste.
if (r_squared > 0.5) {
cat("El R-squared es alto, indicando una buena proporción de varianza explicada por el modelo.\n")
} else {
cat("El R-squared is bajo, lo que sugiere que el modelo puede no capturar gran parte de la varianza en los datos.\n")
}
## El R-squared is bajo, lo que sugiere que el modelo puede no capturar gran parte de la varianza en los datos.
if (adjusted_r_squared > 0.5) {
cat("El Adjusted R-squared es alto, lo cual es favorable al considerar el número de predictores.\n")
} else {
cat("El Adjusted R-squared es bajo, sugiriendo que predictores adicionales pueden mejorar el modelo.\n")
}
## El Adjusted R-squared es bajo, sugiriendo que predictores adicionales pueden mejorar el modelo.