Modelo de regresión basado utilzando datos de las erupciones de volcanes. En primer lugar, cargamos las librerias a utilizar.
suppressMessages(library(caret))
suppressMessages(library(kernlab))
suppressMessages(library(e1071))
suppressMessages(library(sjstats))
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.14
## Current Matrix version is 1.2.12
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
Exploramos los datos “faithful”. ¿Que representan los datos? ¿Cuantas observaciones y variables hay? ¿Que podemos estudiar en basa a estos datos?
#cargamos los datos
data(faithful)
#usar ?faithful para ver la descripción funcional de los datos
#utilizamos este comando para que pueda reproducirse los comandos que utilizan azar
set.seed(1234)
#exploramos las primeras obervaciones
head(faithful)
## eruptions waiting
## 1 3.600 79
## 2 1.800 54
## 3 3.333 74
## 4 2.283 62
## 5 4.533 85
## 6 2.883 55
#revisamos la estructura de los datos
str(faithful)
## 'data.frame': 272 obs. of 2 variables:
## $ eruptions: num 3.6 1.8 3.33 2.28 4.53 ...
## $ waiting : num 79 54 74 62 85 55 88 85 51 85 ...
Creamos las particiones de los datos, utilizando 50% para training y el resto para test.
inTrain<-createDataPartition(y=faithful$waiting,p=0.5,list=FALSE)
training<-faithful[inTrain,]
testing<-faithful[-inTrain,]
dim(training)
## [1] 137 2
dim(testing)
## [1] 135 2
Creamos un modelo de regresión lineal con el tiempo de espera y el tiempo que dura la predicción. Hacemos un diagarma de los tiempos de erupción según el tiempo de espera.
lineal_model<-lm(eruptions~waiting, data=training)
plot(training$waiting, training$eruptions,pch=19, col="blue", xlab="waiting", ylab="Durations")
lines(training$waiting,lineal_model$fitted,lwd=3)
summary(lineal_model)
##
## Call:
## lm(formula = eruptions ~ waiting, data = training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05781 -0.37892 0.04796 0.36478 1.22872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.986918 0.234128 -8.486 3.36e-14 ***
## waiting 0.076647 0.003251 23.579 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5115 on 135 degrees of freedom
## Multiple R-squared: 0.8046, Adjusted R-squared: 0.8032
## F-statistic: 556 on 1 and 135 DF, p-value: < 2.2e-16
Coeficientes del modelo: Indicar cual es son los coeficientes del modelo. ¿Qué significa un coeficiente de 0.076482 para “waiting”? ¿Qué significa un Adjusted R-squares of 0.8056?
lineal_model$coefficients
## (Intercept) waiting
## -1.98691839 0.07664729
summary(lineal_model)
##
## Call:
## lm(formula = eruptions ~ waiting, data = training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05781 -0.37892 0.04796 0.36478 1.22872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.986918 0.234128 -8.486 3.36e-14 ***
## waiting 0.076647 0.003251 23.579 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5115 on 135 degrees of freedom
## Multiple R-squared: 0.8046, Adjusted R-squared: 0.8032
## F-statistic: 556 on 1 and 135 DF, p-value: < 2.2e-16
¿Cual es el tiempo de erupción esperado para una espera de 80 minutos? Podemos calcularlo utilizando los coeficientes o la función “predict”
#calculamos el valor en forma "manual" en base a la recta de regresión
coef(lineal_model)[1]+coef(lineal_model)[2]*80
## (Intercept)
## 4.144865
#Calculamos utilizando la funcion "predict"
predict(lineal_model, newdata=data.frame(waiting=80), interval="prediction")
## fit lwr upr
## 1 4.144865 3.127754 5.161975
Predecimos las erupciones para el set de datos de test
#Predecimos para los datos de testing
pred_eruptions<- predict(lineal_model, newdata=testing,interval="prediction")
head(pred_eruptions)
## fit lwr upr
## 1 4.068217 3.0514624 5.084972
## 5 4.528101 3.5086084 5.547594
## 7 4.758043 3.7366377 5.779448
## 9 1.922093 0.8988047 2.945382
## 10 4.528101 3.5086084 5.547594
## 11 2.152035 1.1309613 3.173109
Calculamos el error RMSE como la suma del cuadrado de las diferencias entre el valor predecido y los valores reales.
#Errores RMSE de la prediccion
sqrt(mean((pred_eruptions-testing$eruptions)^2))
## [1] 0.9631048
Creamos el modelo utilizando Caret y revisamos los resultados - Coeficientes - RMSE - R-squared
modelo_lm <-train(eruptions~waiting, data=training, method="lm")
summary(modelo_lm)
##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05781 -0.37892 0.04796 0.36478 1.22872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.986918 0.234128 -8.486 3.36e-14 ***
## waiting 0.076647 0.003251 23.579 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5115 on 135 degrees of freedom
## Multiple R-squared: 0.8046, Adjusted R-squared: 0.8032
## F-statistic: 556 on 1 and 135 DF, p-value: < 2.2e-16