#==============================================================
# INSTALAR PAQUETES (solo la primera vez)
#==============================================================
paquetes <- c("tidyverse",
"car",
"lmtest",
"sandwich",
"performance",
"GGally"
)
instalar <- paquetes[!(paquetes %in% installed.packages()[,"Package"])]
if(length(instalar)>0){
install.packages(instalar)
}
lapply(paquetes, library, character.only = TRUE)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── 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
## Cargando paquete requerido: carData
##
##
## Adjuntando el paquete: 'car'
##
##
## The following object is masked from 'package:dplyr':
##
## recode
##
##
## The following object is masked from 'package:purrr':
##
## some
##
##
## Cargando paquete requerido: zoo
##
##
## Adjuntando el paquete: 'zoo'
##
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "car" "carData" "lubridate" "forcats" "stringr" "dplyr"
## [7] "purrr" "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[3]]
## [1] "lmtest" "zoo" "car" "carData" "lubridate" "forcats"
## [7] "stringr" "dplyr" "purrr" "readr" "tidyr" "tibble"
## [13] "ggplot2" "tidyverse" "stats" "graphics" "grDevices" "utils"
## [19] "datasets" "methods" "base"
##
## [[4]]
## [1] "sandwich" "lmtest" "zoo" "car" "carData" "lubridate"
## [7] "forcats" "stringr" "dplyr" "purrr" "readr" "tidyr"
## [13] "tibble" "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [19] "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "performance" "sandwich" "lmtest" "zoo" "car"
## [6] "carData" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[6]]
## [1] "GGally" "performance" "sandwich" "lmtest" "zoo"
## [6] "car" "carData" "lubridate" "forcats" "stringr"
## [11] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [16] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [21] "utils" "datasets" "methods" "base"
#==============================================================
# BASE DE DATOS
#==============================================================
datos <- data.frame(
AÑO = c(2018,2019,2020,2021,2022,2023,2024,2025),
PR = c(51.0,
52.6,
48.6,
48.5,
45.9,
44.0,
42.5,
39.5),
AT = c(17.1,
20.7,
23.9,
28.8,
32.2,
41.4,
41.9,
56.9)
)
datos
## AÑO PR AT
## 1 2018 51.0 17.1
## 2 2019 52.6 20.7
## 3 2020 48.6 23.9
## 4 2021 48.5 28.8
## 5 2022 45.9 32.2
## 6 2023 44.0 41.4
## 7 2024 42.5 41.9
## 8 2025 39.5 56.9
#==============================================================
# ESTADÍSTICA DESCRIPTIVA
#==============================================================
summary(datos)
## AÑO PR AT
## Min. :2018 Min. :39.50 Min. :17.10
## 1st Qu.:2020 1st Qu.:43.62 1st Qu.:23.10
## Median :2022 Median :47.20 Median :30.50
## Mean :2022 Mean :46.58 Mean :32.86
## 3rd Qu.:2023 3rd Qu.:49.20 3rd Qu.:41.52
## Max. :2025 Max. :52.60 Max. :56.90
sd(datos$PR)
## [1] 4.430979
sd(datos$AT)
## [1] 13.23426
cor(datos[,2:3])
## PR AT
## PR 1.0000000 -0.9681902
## AT -0.9681902 1.0000000
ggplot(datos,
aes(AÑO,PR))+
geom_line()+
geom_point(size=3)+
theme_bw()
ggplot(datos,
aes(AÑO,AT))+
geom_line()+
geom_point(size=3)+
theme_bw()
#==============================================================
# CORRELACIONES
#==============================================================
cor.test(datos$PR,
datos$AT,
method="pearson")
##
## Pearson's product-moment correlation
##
## data: datos$PR and datos$AT
## t = -9.4781, df = 6, p-value = 7.856e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9944157 -0.8293421
## sample estimates:
## cor
## -0.9681902
cor.test(datos$PR,
datos$AT,
method="spearman")
##
## Spearman's rank correlation rho
##
## data: datos$PR and datos$AT
## S = 166, p-value = 0.0003968
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.9761905
#==============================================================
# MODELO MCO
#==============================================================
modelo <- lm(PR~AT,
data=datos)
summary(modelo)
##
## Call:
## lm(formula = PR ~ AT, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1454 -0.8827 -0.2460 0.6353 2.0824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 57.2277 1.2010 47.649 5.73e-09 ***
## AT -0.3242 0.0342 -9.478 7.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.198 on 6 degrees of freedom
## Multiple R-squared: 0.9374, Adjusted R-squared: 0.927
## F-statistic: 89.83 on 1 and 6 DF, p-value: 7.856e-05
#==============================================================
# ANOVA
#==============================================================
anova(modelo)
## Analysis of Variance Table
##
## Response: PR
## Df Sum Sq Mean Sq F value Pr(>F)
## AT 1 128.831 128.831 89.835 7.856e-05 ***
## Residuals 6 8.604 1.434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(modelo)
## 2.5 % 97.5 %
## (Intercept) 54.2889226 60.1665610
## AT -0.4078478 -0.2404743
#==============================================================
# HC3
#==============================================================
coeftest(modelo,
vcov = vcovHC(modelo,
type="HC3"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 57.22774 1.71633 33.3430 4.843e-08 ***
## AT -0.32416 0.04811 -6.7378 0.0005205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#==============================================================
# SHAPIRO
#==============================================================
shapiro.test(residuals(modelo))
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.89264, p-value = 0.2476
#==============================================================
# BREUSCH PAGAN
#==============================================================
bptest(modelo)
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 0.98794, df = 1, p-value = 0.3202
#==============================================================
# DURBIN WATSON
#==============================================================
dwtest(modelo)
##
## Durbin-Watson test
##
## data: modelo
## DW = 3.1754, p-value = 0.9428
## alternative hypothesis: true autocorrelation is greater than 0
#==============================================================
# DISTANCIA DE COOK
#==============================================================
cooks.distance(modelo)
## 1 2 3 4 5 6
## 0.118435558 0.652693036 0.078554783 0.024050228 0.045230343 0.003583873
## 7 8
## 0.134126139 0.655778047
rstudent(modelo)
## 1 2 3 4 5 6 7
## -0.6638955 3.1723773 -0.7911480 0.5123518 -0.7666794 0.1629424 -1.0780905
## 8
## 0.9319563
#==============================================================
# DIAGNOSTICO
#==============================================================
par(mfrow=c(2,2))
plot(modelo)
par(mfrow=c(1,1))
summary(modelo)$r.squared
## [1] 0.9373923
summary(modelo)$adj.r.squared
## [1] 0.9269577
sigma(modelo)
## [1] 1.197532
datos$Predicho <- predict(modelo)
datos
## AÑO PR AT Predicho
## 1 2018 51.0 17.1 51.68459
## 2 2019 52.6 20.7 50.51761
## 3 2020 48.6 23.9 49.48029
## 4 2021 48.5 28.8 47.89190
## 5 2022 45.9 32.2 46.78976
## 6 2023 44.0 41.4 43.80748
## 7 2024 42.5 41.9 43.64539
## 8 2025 39.5 56.9 38.78298
ggplot(datos,
aes(PR,Predicho))+
geom_point(size=3)+
geom_abline(intercept=0,
slope=1,
color="red")+
theme_bw()
datos$Residuales <- residuals(modelo)
datos
## AÑO PR AT Predicho Residuales
## 1 2018 51.0 17.1 51.68459 -0.6845882
## 2 2019 52.6 20.7 50.51761 2.0823915
## 3 2020 48.6 23.9 49.48029 -0.8802932
## 4 2021 48.5 28.8 47.89190 0.6080958
## 5 2022 45.9 32.2 46.78976 -0.8897567
## 6 2023 44.0 41.4 43.80748 0.1925248
## 7 2024 42.5 41.9 43.64539 -1.1453947
## 8 2025 39.5 56.9 38.78298 0.7170207