Applicazioni pratiche di machine learning/Previsioni sulle produzioni agricole

Caricamento librerie

 library(dplyr)
 library(caret)

Parte 1 : Dati

Il dataset 'agricoltura' scaricabile da qui : https://www.kaggle.com/anjali21/agricultural-production-india contiene dati relativi alle produzioni agricole in India . Esso contiene 246.091 records contenenti ciascuno le seguenti variabili :

• State_Name: Nome dello stato indiano
• District_Name: Nome del distretto all'interno dello stato
• Crop_Year: Anno di produzione del raccolto
• Season: stagione della coltivazione
• Crop: Nome della prodotto coltivato
• Area: Area coltivata
• Production: Produzione del raccolto

Caricamento dei dati:

agricoltura <- read.csv("apy.csv")
agricoltura$Production <- as.numeric(agricoltura$Production)

Parte 2: Esplorazione dati

Si creano 2 datasets relativi alla produzione media di grano e miglio dal 1997 al 2015 e si visualizzano le relative serie storiche:

df1<-agricoltura %>%
  filter(Crop=="Wheat") %>%
  group_by(Crop_Year) %>%
  summarise(Area= mean(Area),Production=mean(Production))
  
df2<-agricoltura %>%
  filter(Crop=="Small millets") %>%
  group_by(Crop_Year) %>%
  summarise(Area= mean(Area),Production=mean(Production))

df1 %>%
  ggplot(aes(x=Crop_Year,y=Production, group=1))+
  geom_line(color="red")+
  geom_text(aes(label=round(Production)), vjust=0, size=2)+
  scale_x_continuous(breaks=df1$Crop_Year)+
  theme(axis.text.x = element_text(angle=45,hjust=1))+
  xlab("Anni")+
  ggtitle("Serie storica del grano in India", subtitle = "dal 1997 al 2015")

df2 %>%
  ggplot(aes(x=Crop_Year,y=Production, group=1))+
  geom_line(color="red")+
  geom_text(aes(label=round(Production)), vjust=0, size=2)+
  scale_x_continuous(breaks=df1$Crop_Year)+
  theme(axis.text.x = element_text(angle=45,hjust=1))+
  xlab("Anni")+
  ggtitle("Serie storica del miglio in India", subtitle = "dal 1997 al 2015")

Parte 3 : Modellizzazione e previsione

Si utizza una regressione polinomiale di grado 5 per prevedere la produzione di grano negli anni successivi al 2015, ottenendo un R quadro del modello pari all'81,3%

n<-nrow(df1)
q <- 1:n
df <- data.frame(production=df1$Production,x=poly(q,5, raw=TRUE))
model <- lm(production ~ .,df)
summary(model)

Call:
lm(formula = production ~ ., data = df)

Residuals:
  Min     1Q Median     3Q    Max 
-983.8 -418.0   71.3  357.0 1103.3 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.714e+04  1.559e+03  17.415 2.16e-10 ***
x.1         -2.101e+03  1.428e+03  -1.471  0.16500    
x.2          8.809e+02  4.152e+02   2.122  0.05365 .  
x.3         -1.397e+02  5.110e+01  -2.735  0.01702 *  
x.4          9.232e+00  2.787e+00   3.312  0.00561 ** 
x.5         -2.137e-01  5.551e-02  -3.851  0.00200 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 662.9 on 13 degrees of freedom
Multiple R-squared:  0.8649,	Adjusted R-squared:  0.813 
F-statistic: 16.65 on 5 and 13 DF,  p-value: 2.999e-05

Dopo 2 anni secondo il modello regressivo la produzione di grano si dovrebbe azzerare :

df_pred <-data.frame(x=poly(20:21,5, raw=TRUE))
predict.lm (model,newdata=df_pred, interval="confidence", level=0.80)

         fit       lwr       upr
1 12645.2107 10540.908 14749.513
2  -204.5818 -4776.246  4367.082

plot(1:n,df1$Production, main = "Serie storica del grano con polinomio di 5° grado", col="blue", type="l" , xlab = "Anni",ylab =  "Produzione")
  points(1:n,predict(model,df), col="red",type="b", lty=2)
  legend(1,23000,legend = c("serie storica","previsione"),col = c("blue","red"), lty = 1:2, cex = 0.8)

Si utizza una regressione polinomiale di grado 6 per prevedere la produzione di miglio negli anni successivi al 2015, ottenendo un R quadro del modello pari al 90,37%

n<-nrow(df2)
q <- 1:n
df <- data.frame(production=df2$Production,x=poly(q,6, raw=TRUE))
model <- lm(production ~ .,df)
summary(model)

Call:
lm(formula = production ~ ., data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-3022.94  -525.28    -8.86   991.11  2183.66  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.202e+04  5.175e+03   6.188 4.67e-05 ***
x.1         -1.380e+04  6.181e+03  -2.233 0.045393 *  
x.2          6.588e+03  2.442e+03   2.698 0.019375 *  
x.3         -1.364e+03  4.354e+02  -3.133 0.008643 ** 
x.4          1.372e+02  3.861e+01   3.553 0.003974 ** 
x.5         -6.563e+00  1.660e+00  -3.952 0.001919 ** 
x.6          1.196e-01  2.760e-02   4.333 0.000974 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1491 on 12 degrees of freedom
Multiple R-squared:  0.9358,	Adjusted R-squared:  0.9037 
F-statistic: 29.14 on 6 and 12 DF,  p-value: 1.754e-06

la produzione di miglio secondo il modello previsionale tenderà negli anni successivi a crescere nel seguente modo con un intervallo di confidenza all'80% :

df_pred <-data.frame(x=poly(20:30,6, raw=TRUE))
predict.lm (model,newdata=df_pred, interval="confidence", level=0.80)

          fit        lwr        upr
1    77545.92   70527.07   84564.77
2   143944.02  124892.30  162995.75
3   260084.67  218097.79  302071.55
4   449274.05  367753.40  530794.69
5   741632.36  596404.25  886860.47
6  1175243.26  932368.21 1418118.32
7  1797389.41 1410612.72 2184166.10
8  2665874.07 2073687.17 3258060.98
9  3850428.88 2972713.36 4728144.41
10 5434207.62 4168434.96 6699980.29
11 7515366.14 5732325.90 9298406.38

plot(1:n,df2$Production, main = "Serie storica del miglio in India", col="blue", type="l" , xlab = "Anni",ylab =  "Produzione")
  points(1:n,predict(model,df), col="red",type="b", lty=2)
  legend(1,35000,legend = c("serie storica","previsione"),col = c("blue","red"), lty = 1:2, cex = 0.8)