Una cooperativa agrícula quiere estudiar los factores que influiyen en la producción de maiz de hectárea.
library(tidyverse)
set.seed(1234)
datos <- data.frame(
Fertilizacion = round(runif(100, 20, 100)), #kg/ha
Riego = round(runif(100, 50, 200)), #mm de agua
Plagas = round(runif(100, 0, 10)), #indice de plagas
error = rnorm(100, mean = 0, sd = 5)
) %>%
mutate(Producción = 50 + (2 * Fertilizacion) + (1.5 * Riego) - (3 * Plagas) + error)
head(datos)
## Fertilizacion Riego Plagas error Producción
## 1 29 55 7 -1.8861882 167.6138
## 2 70 135 5 0.4880973 377.9881
## 3 69 92 3 8.1937232 325.1937
## 4 70 81 8 -4.3779624 283.1220
## 5 89 70 5 0.6088000 318.6088
## 6 71 99 7 6.8106533 326.3107
Utilizamos la herramienta de set.seed() para generar números aleatorios y generamos un conjunto de de daros con 100 observaciones para las variables.
# Ajustar modelo de regresión lineal múltiple
modelo <- lm(Producción ~ Fertilizacion + Riego + Plagas, data = datos)
# Resumen del modelo
summary(modelo)
##
## Call:
## lm(formula = Producción ~ Fertilizacion + Riego + Plagas, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.2699 -2.6389 0.1072 2.5374 14.5880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.93009 2.37763 20.58 <2e-16 ***
## Fertilizacion 2.00345 0.02517 79.60 <2e-16 ***
## Riego 1.50934 0.01261 119.66 <2e-16 ***
## Plagas -2.92918 0.20683 -14.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.503 on 96 degrees of freedom
## Multiple R-squared: 0.9952, Adjusted R-squared: 0.995
## F-statistic: 6570 on 3 and 96 DF, p-value: < 2.2e-16
El valor de p es sumamente bajo significando que la fertilización, el riego y las plagas tienen un efecto significativo en la producción. Además el valor de R^2 es 0.995, significando que el modelo es sumamente preciso, demostrando la variabilidad entre la producción y las variables simuladas.
Lote A: fertilización = 40, riego = 100, plagas = 2, Lote B: fertilización = 80, riego = 150, plagas = 5
nuevos <- data.frame(
Fertilizacion=c(40,80),
Riego=c(100,150),
Plagas=c(2,5)
)
#Predicción puntual
round(predict(modelo, nuevos))
## 1 2
## 274 421
#Intervalo de confianza (IC 92%)
round(predict(modelo, nuevos, interval = "confidence", level = 0.92))
## fit lwr upr
## 1 274 272 276
## 2 421 419 423
#Intervalo de predicción (IP 98%)
round(predict(modelo, nuevos, interval = "prediction", level = 0.98))
## fit lwr upr
## 1 274 261 287
## 2 421 408 434
Se evaluaron dos casos diferentes, Lote A y Lote B.
Lote A Lote A proyecta un rendimiento de 274kg/ha con un 92% de confianza entre 274 kg/ha y 276 kg/ha y con un intervalo de predicción de 98% de 261 kg/ha y 287 kg/ha.
Lote B proyecta un rendimiento de 421 kg/ha con un 92% de confianza entre 419 kg/ha y 423 kg/ha y con un intervalo de predicción de 98% de 408 kg/ha a 434 kg/ha.
library(tidyverse)
library(plotly)
# Escenarios de simulación
simulacion1 <- data.frame(
Plagas = seq(0, 10, by = 2),
Fertilizacion = 60,
Riego = 120)
# Predicciones con IC
pred_conf <- data.frame(predict(modelo, newdata = simulacion1, interval = "confidence",level = 0.92))
colnames(pred_conf) <- c("fit", "lwr_conf", "upr_conf")
# Predicciones con IP
pred_pred <- data.frame(predict(modelo, newdata = simulacion1, interval = "prediction",level = 0.98))
colnames(pred_pred) <- c("fit_pred", "lwr_pred", "upr_pred")
simulacion2 <- bind_cols(simulacion1, pred_conf, pred_pred)
# Gráfico interactivo con plotly
p <- ggplot() +
# Observaciones reales
geom_point(data = datos, aes(x = Plagas, y = Producción),
color = "red", size = 1.2, alpha = 0.7) +
# Línea de predicción puntual
geom_line(data = simulacion2, aes(x = Plagas, y = fit),
color = "green4", size = 0.5) +
# Banda de IC
geom_ribbon(data = simulacion2, aes(x = Plagas, ymin = lwr_conf, ymax = upr_conf),
fill = "purple", alpha = 0.3) +
# Banda de IP
geom_ribbon(data = simulacion2, aes(x = Plagas, ymin = lwr_pred, ymax = upr_pred),
fill = "pink", alpha = 0.15) +
labs(title = "",
x = "Indice de Plagas",
y = "Producción") +
theme_minimal()
ggplotly(p)
Aqui podemos ver las 100 observaciones muestran cómo disminuye la producción a medida que aumentan las plagas, manteniendo fijos la fertilización en 60 y el riego en 120.