Generar el modelo

Se utilizará el método mixto iniciando el modelo con todas las variables como predictores y realizando la selección de los mejores predictores con la medición Akaike (AIC).

modelo <- lm(ProduccionCarne~
               Precipitaciones + Temperaturas , data = datos )
summary(modelo)

## 
## Call:
## lm(formula = ProduccionCarne ~ Precipitaciones + Temperaturas, 
##     data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65738 -35814 -15579  23269 187877 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)      54426.94  116745.77   0.466    0.645
## Precipitaciones     20.38      24.02   0.849    0.403
## Temperaturas      -308.05    4054.68  -0.076    0.940
## 
## Residual standard error: 61660 on 29 degrees of freedom
## Multiple R-squared:  0.02534,    Adjusted R-squared:  -0.04188 
## F-statistic: 0.377 on 2 and 29 DF,  p-value: 0.6892

El modelo con todas las variables introducidas como predictores tiene un R^2 alta (0.02534), es capaz de explicar el 3% de la variabilidad observada en las precipitaciones. El p-value del modelo es significativo (0.6892) por lo que se puede aceptar que el modelo no es por azar, al menos uno de los coeficientes parciales de regresión es distinto de 0.

Selección de los mejores predictores

step(object = modelo, direction = "both", trace = 1)

## Start:  AIC=708.73
## ProduccionCarne ~ Precipitaciones + Temperaturas
## 
##                   Df  Sum of Sq        RSS    AIC
## - Temperaturas     1   21942743 1.1027e+11 706.73
## - Precipitaciones  1 2737643509 1.1298e+11 707.51
## <none>                          1.1024e+11 708.73
## 
## Step:  AIC=706.73
## ProduccionCarne ~ Precipitaciones
## 
##                   Df  Sum of Sq        RSS    AIC
## - Precipitaciones  1 2844468357 1.1311e+11 705.55
## <none>                          1.1027e+11 706.73
## + Temperaturas     1   21942743 1.1024e+11 708.73
## 
## Step:  AIC=705.55
## ProduccionCarne ~ 1
## 
##                   Df  Sum of Sq        RSS    AIC
## <none>                          1.1311e+11 705.55
## + Precipitaciones  1 2844468357 1.1027e+11 706.73
## + Temperaturas     1  128767591 1.1298e+11 707.51

## 
## Call:
## lm(formula = ProduccionCarne ~ 1, data = datos)
## 
## Coefficients:
## (Intercept)  
##       63364

modelo <- lm(ProduccionCarne~
               Precipitaciones + Temperaturas , data = datos )
summary(modelo)

## 
## Call:
## lm(formula = ProduccionCarne ~ Precipitaciones + Temperaturas, 
##     data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65738 -35814 -15579  23269 187877 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)      54426.94  116745.77   0.466    0.645
## Precipitaciones     20.38      24.02   0.849    0.403
## Temperaturas      -308.05    4054.68  -0.076    0.940
## 
## Residual standard error: 61660 on 29 degrees of freedom
## Multiple R-squared:  0.02534,    Adjusted R-squared:  -0.04188 
## F-statistic: 0.377 on 2 and 29 DF,  p-value: 0.6892

confint(lm(ProduccionCarne~
               Precipitaciones + Temperaturas , data = datos ))

##                         2.5 %       97.5 %
## (Intercept)     -184344.95758 293198.84206
## Precipitaciones     -28.73707     69.49626
## Temperaturas      -8600.80102   7984.69954

Validación de condiciones para la regresión múltiple lineal

library(gridExtra)

## Loading required package: grid

ggplot(data = datos, aes(Precipitaciones, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()

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

ggplot(data = datos, aes(Temperaturas, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()

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

Download bovinoNacional2019.xlsx

Distribución normal de los residuos:

qqnorm(modelo$residuals)
qqline(modelo$residuals)

Download bovinoNacional2019.xlsx

Shapiro test:

shapiro.test(modelo$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.81051, p-value = 6.464e-05

Tanto el análisis gráfico como es test de hipótesis (Shapiro test) confirman la normalidad.

ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()

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

Download bovinoNacional2019.xlsx

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.0.5

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.0.5

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(modelo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 0.8354, df = 2, p-value = 0.6586

Matriz de correlación entre predictores

library(corrplot)

## Warning: package 'corrplot' was built under R version 4.0.5

## corrplot 0.88 loaded

corrplot(cor(dplyr::select(datos,Precipitaciones,Temperaturas)),
         method = "number", tl.col = "black")

Download bovinoNacional2019.xlsx

Análisis de Inflación de Varianza (VIF):

library(car)

## Warning: package 'car' was built under R version 4.0.5

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

vif(modelo)

## Precipitaciones    Temperaturas 
##        1.096488        1.096488

Autocorrelación:

library(car)
dwt(modelo, alternative = "two.sided")

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.2321819      2.461288   0.146
##  Alternative hypothesis: rho != 0

Identificación de posibles valores atípicos o influyentes

library(dplyr)
datos$studentized_residual <- rstudent(modelo)
ggplot(data = datos, aes(x = predict(modelo), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
scale_color_identity() +
labs(title = "Distribución",
     x = "predicción modelo") + 
theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Download bovinoNacional2019.xlsx

which(abs(datos$studentized_residual) > 3)

## [1] 14 30

summary(influence.measures(modelo))

## Potentially influential observations of
##   lm(formula = ProduccionCarne ~ Precipitaciones + Temperaturas,      data = datos) :
## 
##    dfb.1_ dfb.Prcp dfb.Tmpr dffit   cov.r   cook.d hat    
## 5   0.04   0.18    -0.06     0.19    1.40_*  0.01   0.22  
## 14 -0.05  -0.10     0.12     0.63    0.40_*  0.10   0.03  
## 15 -0.17  -0.04     0.16    -0.18    1.34_*  0.01   0.18  
## 27  0.05  -0.26    -0.01    -0.29    1.66_*  0.03   0.34_*
## 30  0.37   1.05_*  -0.48     1.29_*  0.33_*  0.37   0.10

La visualización gráfica de las influencias se obtiene del siguiente modo:

influencePlot(modelo)

##       StudRes        Hat      CookD
## 5   0.3649458 0.22001761 0.01290884
## 14  3.4332368 0.03279691 0.09710847
## 27 -0.4015391 0.34096410 0.02863389
## 30  3.9236591 0.09750239 0.37049963

Download bovinoNacional2019.xlsx

\[ProduccionCarne = 20.38 Precipitaciones -308.05 Temperaturas \] Es capaz de explicar el 2.53% de la variabilidad observada en las precipitaciones (R2: 0.02534, R2-Adjusted: -0.04188). El test F muestra que es significativo (p-value: 0.6892). Se satisfacen todas las condiciones para este tipo de regresión múltiple.Lo que significa que las observaciones podrían estar influyendo de forma notable en el modelo.