Agronomic Growth2

library(tidyverse)
library(modelr)
library(dplyr)

library(readxl)
Agronomic_Yields <- read_excel("~/Agronomic Yields.xlsx")
View(Agronomic_Yields)

Agronomic_Yields_Model2 <- lm(formula = yield_kg_ha ~ treatment, data = Agronomic_Yields)

summary(Agronomic_Yields_Model2)

Call:
lm(formula = yield_kg_ha ~ treatment, data = Agronomic_Yields)

Residuals:
    Min      1Q  Median      3Q     Max 
-5706.7 -2582.7  -792.5  2276.2  7814.4 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4781.34     205.95   23.22   <2e-16 ***
treatment     354.04      36.61    9.67   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3244 on 1178 degrees of freedom
  (16 observations deleted due to missingness)
Multiple R-squared:  0.07355,   Adjusted R-squared:  0.07276 
F-statistic: 93.52 on 1 and 1178 DF,  p-value: < 2.2e-16

Agronomic_Yields_Model4 <- lm(formula = yield_kg_ha ~ fertilizer_rate_kg_ha, data = Agronomic_Yields)

summary(Agronomic_Yields_Model4)

Call:
lm(formula = yield_kg_ha ~ fertilizer_rate_kg_ha, data = Agronomic_Yields)

Residuals:
    Min      1Q  Median      3Q     Max 
-7598.2 -2034.6  -580.7  2009.8  7872.8 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           4437.6658   140.8019   31.52   <2e-16 ***
fertilizer_rate_kg_ha   20.0739     0.9933   20.21   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2951 on 1035 degrees of freedom
  (159 observations deleted due to missingness)
Multiple R-squared:  0.2829,    Adjusted R-squared:  0.2823 
F-statistic: 408.4 on 1 and 1035 DF,  p-value: < 2.2e-16

Agronomic_Yields_Model5 <- lm(formula = yield_kg_ha ~ fertilizer_rate_kg_ha + treatment, data = Agronomic_Yields)

summary(Agronomic_Yields_Model5)

Call:
lm(formula = yield_kg_ha ~ fertilizer_rate_kg_ha + treatment, 
    data = Agronomic_Yields)

Residuals:
   Min     1Q Median     3Q    Max 
 -7583  -2319   -363   1917   7171 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5772.25     197.39  29.242   <2e-16 ***
fertilizer_rate_kg_ha    30.94       1.51  20.489   <2e-16 ***
treatment              -501.27      53.97  -9.288   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2836 on 1034 degrees of freedom
  (159 observations deleted due to missingness)
Multiple R-squared:  0.3382,    Adjusted R-squared:  0.3369 
F-statistic: 264.2 on 2 and 1034 DF,  p-value: < 2.2e-16

Agronomic_Yields_Model3 <- lm(formula = yield_kg_ha ~ fertilizer_rate_kg_ha + crop, data = Agronomic_Yields)

summary(Agronomic_Yields_Model3)

Call:
lm(formula = yield_kg_ha ~ fertilizer_rate_kg_ha + crop, data = Agronomic_Yields)

Residuals:
    Min      1Q  Median      3Q     Max 
-7880.9 -1354.0    84.9  1157.5  5767.9 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  3386.4662   161.0970  21.021  < 2e-16 ***
fertilizer_rate_kg_ha          12.9836     0.8993  14.438  < 2e-16 ***
cropTriticum aestivum L. (*) -916.5881   238.4176  -3.844 0.000128 ***
cropZea mays L. (*)          3397.2160   212.2105  16.009  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2337 on 1033 degrees of freedom
  (159 observations deleted due to missingness)
Multiple R-squared:  0.5513,    Adjusted R-squared:   0.55 
F-statistic:   423 on 3 and 1033 DF,  p-value: < 2.2e-16

ggplot(data = Agronomic_Yields,
       mapping = aes(x = as.factor(treatment), y = yield_kg_ha)) +
stat_summary(fun = "mean",
             geom = "bar",
             fill = "grey",
             color = "black") +
  stat_summary(fun.data = "mean_se", 
               geom = "errorbar",
               width = 0.2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank()) +
labs(x = "Treatment",
     y = "Yield (kg/ha)")

ggplot(data = Agronomic_Yields,
       mapping = aes(x = fertilizer_rate_kg_ha, y = yield_kg_ha, color = factor(treatment))) +
  stat_summary(fun = "mean", geom = "point", size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank()) +
labs(x = "Fertilizer Rate (kg/ha)",
     y = "Yield (kg/ha)",
     color = "Treatment")

ggplot(data = Agronomic_Yields,
       mapping = aes(x = fertilizer_rate_kg_ha, y = yield_kg_ha)) +
  stat_summary(fun = "mean", geom = "point", size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank()) +
labs(x = "Fertilizer Rate (kg/ha)",
     y = "Yield (kg/ha)")

ggplot(data = Agronomic_Yields, 
       mapping = aes(x = fertilizer_rate_kg_ha, y = yield_kg_ha, color = crop)) +
  stat_summary(fun = "mean", geom = "point", size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank()) +
labs(x = "Fertilizer Rate (kg/ha)",
     y = "Yield (kg/ha)",
     color = "Crop Type")