Introduction

In this project, prepared as part of the course IE 2152 – Statistics for Industrial Engineers, a statistical analysis has been conducted using the agricultural production data provided by the Food and Agriculture Organization (FAO). According to the student number rule, I selected Argentina as the target country based on the last digit of my student ID (0 = A).

From Argentina’s agricultural data, Apple has been chosen as the product of interest, and the analysis focuses on its Gross Production Value (constant 2014-2016 thousand I$) between 1961 and 2023.

Throughout the project, the data set will be examined under the following four main sections:

  1. Descriptive Statistics
  2. Inferential Statistics
  3. Exploratory Data Analysis
  4. Trend-Based Forecasting

Each section includes appropriate statistical techniques, visualizations, and interpretations to provide a comprehensive analysis of the selected data.

Section 0-Determining the Elements to be Used

# We need to read the CSV file from my local path
apple_raw <- read.csv("~/Downloads/FAOSTAT_data_en_6-13-2025.csv")
# Viewing column names
colnames(apple_raw)
##  [1] "Domain.Code"      "Domain"           "Area.Code..M49."  "Area"            
##  [5] "Element.Code"     "Element"          "Item.Code..CPC."  "Item"            
##  [9] "Year.Code"        "Year"             "Unit"             "Value"           
## [13] "Flag"             "Flag.Description"

Monitoring Apple’s Years:

years <- apple_raw$Year ;years
##  [1] 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975
## [16] 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990
## [31] 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
## [46] 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
## [61] 2021 2022 2023

Monitoring Apple’s Values:

values <- apple_raw$Value ; values
##  [1] 195881 187385 223729 175113 256769 195409 243553 221652 205981 210041
## [11] 200081 241806 110118 370993 286977 272109 387041 382321 458785 452177
## [21] 428577 379489 385625 411585 435374 280322 507213 443681 400163 460485
## [31] 503862 497254 448954 474975 540914 575370 527551 487823 526754 393329
## [41] 674397 546025 617123 595874 569333 519202 472001 448401 421969 339841
## [51] 542802 495602 457841 335618 327227 216605 216605 240946 258584 277065
## [61] 242609 199819 249570

Section 1-Descriptive Statistics

1.1-Construct frequency distribution, cumulative frequency plot and histogram for the selected data set.

Constructing Frequency Distribution

# Define class width and breaks
breaks <- seq(from = min(values), to = max(values) + 70550, by = 70550)

# Cut the values into class intervals
value.cut <- cut(values, breaks = breaks, right = FALSE)

# Frequency table
value.freq <- table(value.cut)

# Relative frequency
value.relfreq <- value.freq / length(values)

# Cumulative frequency
value.cumfreq <- cumsum(value.freq)

# Cumulative relative frequency
value.cumrelfreq <- value.cumfreq / length(values)

To construct the full frequency table, I first calculated the range of the data by subtracting the minimum value from the maximum. Based on the sample size (n = 63), I determined the optimal number of intervals (k = 8) using the rule \(2^k \geq n\). Then, the class width was obtained by dividing the range by the number of intervals. Using this class width, the data was grouped into bins, and absolute, relative, cumulative, and cumulative relative frequencies were calculated accordingly.

# Display all frequencies in a single table
options(digits = 3)
cbind(value.freq, value.relfreq, value.cumfreq, value.cumrelfreq)
##                     value.freq value.relfreq value.cumfreq value.cumrelfreq
## [1.1e+05,1.81e+05)           2        0.0317             2           0.0317
## [1.81e+05,2.51e+05)         16        0.2540            18           0.2857
## [2.51e+05,3.22e+05)          6        0.0952            24           0.3810
## [3.22e+05,3.92e+05)          8        0.1270            32           0.5079
## [3.92e+05,4.63e+05)         13        0.2063            45           0.7143
## [4.63e+05,5.33e+05)         10        0.1587            55           0.8730
## [5.33e+05,6.04e+05)          6        0.0952            61           0.9683
## [6.04e+05,6.75e+05)          2        0.0317            63           1.0000

Cumulative Frequency Plot

cum_freq <- cumsum(value.freq)
plot(cum_freq, type = "o", col = "red", xlab = "Interval Index", ylab = "Cumulative Frequency",
     main = "Cumulative Frequency Plot")

Histogram

hist(values, breaks = 8, col = "lightgreen", border = "black",
     main = "Histogram of Apple Gross Production Value",
     xlab = "Production Value (thousand I$)", ylab = "Frequency")

1.2-Evaluate central tendency and dispersion measures. Find the mean, standard deviation, variance, median, mode, IQR, and range for the selected dataset.

Mean:

mean(values)
## [1] 375433

The mean value of apple production is approximately 375433 thousand I$.

Standart Deviation:

sd(values)
## [1] 136537

The standard deviation is about 136537, indicating how much values vary around the mean.

Variance:

var(values)
## [1] 1.86e+10

The variance is 1.86e+10, representing the overall dispersion in the dataset.

Median:

median(values)
## [1] 387041

The median value is 387041, showing the central point of the distribution.

Mode:

sort(values)
##  [1] 110118 175113 187385 195409 195881 199819 200081 205981 210041 216605
## [11] 216605 221652 223729 240946 241806 242609 243553 249570 256769 258584
## [21] 272109 277065 280322 286977 327227 335618 339841 370993 379489 382321
## [31] 385625 387041 393329 400163 411585 421969 428577 435374 443681 448401
## [41] 448954 452177 457841 458785 460485 472001 474975 487823 495602 497254
## [51] 503862 507213 519202 526754 527551 540914 542802 546025 569333 575370
## [61] 595874 617123 674397

The mode of the dataset is 216,605, which appears twice. All other values occur only once, indicating that the dataset is nearly uniform in terms of frequency.

Interquartile Range (IQR):

IQR(values)
## [1] 238318

The interquartile range (IQR) is 238318, reflecting the spread of the middle 50% of the data.

Range:

range(values)          # to see both min and max
## [1] 110118 674397
diff(range(values))    # just the numeric range
## [1] 564279

The range of the dataset is 564279, showing the difference between the maximum and minimum production values.

1.3-Find 5 summary statistics of box and whisker plot, respectively

Five-Number Summary

summary(values)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  110118  243081  387041  375433  481399  674397

The five-number summary is as follows:

  • Minimum: 110118

  • First Quartile (Q1): 243081

  • Median: 387041

  • Third Quartile (Q3): 481399

  • Maximum: 674397

These five values summarize the spread and central tendency of the dataset, and form the basis of the boxplot visualization.

Box and Whisker Plot:

boxplot(values,
        main = "Boxplot of Apple Gross Production Value",
        ylab = "Production Value (thousand I$)",
        col = "yellow")

Section 2-Inferential Statistics

2.1 – Hypothesis Testing on a Random Sample (n = 30)

Step 1: Drawing a Random Sample (n = 30)

set.seed(123)  # For reproducibility
sample30 <- sample(values, size = 30)

A random sample of size 30 was drawn from the population using the sample() function. A fixed seed was set to ensure that results remain reproducible.

Step 2: Defining Population Parameters

pop_mean <- mean(values)
pop_var <- var(values)

The population mean and variance were calculated based on the full dataset, and will serve as the null hypothesis reference values.

Step 3: Hypothesis Definition Null Hypothesis (H₀): μ = population mean

Alternative Hypothesis (H₁): μ ≠ population mean (two-tailed)

A two-sided t-test was selected since we do not assume a specific direction of change in the sample mean compared to the population mean.

Step 4: Justification for Using the t-test

Although the population standard deviation was calculated in Part 1, it is not considered “known” in the statistical sense. Therefore, a one-sample t-test is used instead of a z-test.

Step 5: Conducting t-tests for α = 0.01, 0.05, 0.10

Step 5.1 – t-test for α = 0.01 (conf.level = 0.99)

# Define problem data
n <- length(sample30) ;n                       # sample size
## [1] 30
xbar <- mean(sample30) ; xbar                      # sample mean
## [1] 397719
sd <- sd(sample30)   ;sd                        # sample standard deviation
## [1] 128825
mu0 <- pop_mean ;mu0                              # hypothesized population mean
## [1] 375433
alpha <- 0.01                                # significance level for part 5.1
# Calculate the critical t-value for α = 0.01 (two-tailed)
tc <- qt(alpha / 2, df = n - 1, lower.tail = FALSE)
c(-tc, tc)                                   # critical region bounds
## [1] -2.76  2.76
# Calculate the test statistic
t0 <- (xbar - mu0) / (sd / sqrt(n)); t0
## [1] 0.948
plotmean.hypothesis <- function(t0, tc, df, type = "t", alt = "t") {
  x <- seq(-4, 4, length = 1000)
  y <- dt(x, df)
  
  plot(x, y, type = "l", lwd = 2,
       main = "t-distribution with Critical Regions",
       xlab = "t", ylab = "Density")
  
  abline(v = t0, col = "blue", lwd = 2)
  abline(v = c(-tc, tc), col = "red", lty = 2)

  legend("topright", legend = c("Test Statistic", "Critical Values"),
         col = c("blue", "red"), lty = c(1, 2), lwd = 2)
  
  # Shade critical regions
  polygon(c(-4, seq(-tc, -4, length = 100), -tc), 
          c(0, dt(seq(-tc, -4, length = 100), df), 0), col = rgb(1, 0, 0, 0.2), border = NA)
  
  polygon(c(tc, seq(tc, 4, length = 100), 4), 
          c(0, dt(seq(tc, 4, length = 100), df), 0), col = rgb(1, 0, 0, 0.2), border = NA)
}

# Plotting the t-distribution with test statistic

plotmean.hypothesis(t0, tc, n - 1, type = "t", alt = "t")

# Calculate the p-value
pval <- 2 * pt(-abs(t0), df = n - 1); pval
## [1] 0.351
test.pval <- function(pval, mu0, alpha, theta = "mu") {
  cat("Significance level (alpha):", alpha, "\n")
  cat("P-value:", round(pval, 4), "\n")
  
  if (pval < alpha) {
    cat("Decision: Reject H0. There is significant evidence that", theta, "≠", mu0, "\n")
  } else {
    cat("Decision: Fail to reject H0. There is not enough evidence that", theta, "≠", mu0, "\n")
  }
}

# Decision based on p-value
test.pval(pval, mu0, alpha, theta = "mu")
## Significance level (alpha): 0.01 
## P-value: 0.351 
## Decision: Fail to reject H0. There is not enough evidence that mu ≠ 375433

At a significance level of 0.01, the sample mean does not significantly differ from the population mean. The test statistic falls well within the acceptance region, and the p-value (0.351) confirms that we fail to reject the null hypothesis. Therefore, there is insufficient evidence to suggest a meaningful difference at the 1% level of significance.

Step 5.2– t-test for α = 0.05 (conf.level = 0.95)

alpha <- 0.05                                # significance level for part 5.2
# Calculate the critical t-value for α = 0.01 (two-tailed)
tc <- qt(alpha / 2, df = n - 1, lower.tail = FALSE)
c(-tc, tc)                                   # critical region bounds
## [1] -2.05  2.05
# Calculate the test statistic
t0 <- (xbar - mu0) / (sd / sqrt(n)); t0
## [1] 0.948
#same function again

# Plotting the t-distribution with test statistic

plotmean.hypothesis(t0, tc, n - 1, type = "t", alt = "t")

# Calculate the p-value
pval <- 2 * pt(-abs(t0), df = n - 1); pval
## [1] 0.351
#same p value function again
# Decision based on p-value
test.pval(pval, mu0, alpha, theta = "mu")
## Significance level (alpha): 0.05 
## P-value: 0.351 
## Decision: Fail to reject H0. There is not enough evidence that mu ≠ 375433

At a significance level of 0.05, the test still fails to provide strong enough evidence to reject the null hypothesis. Although the confidence level is lower than in the previous test, the test statistic remains within the non-critical region, and the p-value is greater than 0.05. Thus, we conclude that the sample mean is not significantly different from the population mean at the 5% level.

Step 5.3– t-test for α = 0.10 (conf.level = 0.90)

alpha <- 0.10                                # significance level for part 5.3
# Calculate the critical t-value for α = 0.01 (two-tailed)
tc <- qt(alpha / 2, df = n - 1, lower.tail = FALSE)
c(-tc, tc)                                   # critical region bounds
## [1] -1.7  1.7
# Calculate the test statistic
t0 <- (xbar - mu0) / (sd / sqrt(n)); t0
## [1] 0.948
#same function again

# Plotting the t-distribution with test statistic

plotmean.hypothesis(t0, tc, n - 1, type = "t", alt = "t")

# Calculate the p-value
pval <- 2 * pt(-abs(t0), df = n - 1); pval
## [1] 0.351
#same p value function again
# Decision based on p-value
test.pval(pval, mu0, alpha, theta = "mu")
## Significance level (alpha): 0.1 
## P-value: 0.351 
## Decision: Fail to reject H0. There is not enough evidence that mu ≠ 375433

At a significance level of 0.10, we again fail to reject the null hypothesis. While the acceptance region is wider due to a higher alpha, the test statistic is still not extreme enough to indicate a significant difference. The p-value remains above the threshold, supporting the decision that there is no statistically significant difference between the sample and population means at the 10% level.

2.2 – Confidence Intervals (CI)

Step 1 – CI for α = 0.01

ci_99 <- t.test(sample30, conf.level = 0.99)
ci_99$conf.int
## [1] 332888 462549
## attr(,"conf.level")
## [1] 0.99

Step 2 – CI for α = 0.05

ci_95 <- t.test(sample30, conf.level = 0.95)
ci_95$conf.int
## [1] 349615 445823
## attr(,"conf.level")
## [1] 0.95

Step 3 – CI for α = 0.1

ci_90 <- t.test(sample30, conf.level = 0.90)
ci_90$conf.int
## [1] 357755 437682
## attr(,"conf.level")
## [1] 0.9

Step 4 – Prepare Data for Plot

library(ggplot2)


ci_df <- data.frame(
  Confidence_Level = c("99%", "95%", "90%"),
  Lower = c(ci_99$conf.int[1], ci_95$conf.int[1], ci_90$conf.int[1]),
  Upper = c(ci_99$conf.int[2], ci_95$conf.int[2], ci_90$conf.int[2]),
  Mean = mean(sample30)
)

Step 5 – Plot the Confidence Intervals

ggplot(ci_df, aes(x = Confidence_Level, y = Mean)) +
  geom_point(size = 3, color = "blue") +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1, color = "darkred") +
  labs(title = "Confidence Intervals for Sample Mean",
       x = "Confidence Level",
       y = "Production Value (thousand I$)") +
  theme_minimal()

The plot above shows the confidence intervals for the sample mean at three different confidence levels. As expected, the 99% interval is the widest, followed by 95% and 90%. All intervals contain the population mean, further supporting the results of the hypothesis tests.

Step 6 – Create and Print CI Table

ci_df_rounded <- data.frame(
  `Confidence Level` = c("99%", "95%", "90%"),
  `Lower Bound` = round(c(ci_99$conf.int[1], ci_95$conf.int[1], ci_90$conf.int[1]), 2),
  `Mean` = round(mean(sample30), 2),
  `Upper Bound` = round(c(ci_99$conf.int[2], ci_95$conf.int[2], ci_90$conf.int[2]), 2)
)

knitr::kable(ci_df_rounded, caption = "Confidence Intervals for Sample Mean at Different Levels")
Confidence Intervals for Sample Mean at Different Levels
Confidence.Level Lower.Bound Mean Upper.Bound
99% 332888 397719 462549
95% 349615 397719 445823
90% 357755 397719 437682

2.3 – Type I Error

Type I Error – Definition and Interpretation

Type I Error occurs when the null hypothesis is rejected even though it is actually true. In other words, it represents a false positive conclusion. The probability of making a Type I Error is denoted by α (alpha), which is also called the significance level.

The characteristics of Type I Error include: - It is set by the researcher before conducting the test. - It defines the critical region: if the test statistic falls within this region, H₀ is rejected. - The smaller the α, the lower the probability of rejecting a true null hypothesis.

In this project, the hypothesis test was performed using three different significance levels. The corresponding Type I Error probabilities are:

  • α = 0.01 → 1% chance of Type I Error
  • α = 0.05 → 5% chance of Type I Error
  • α = 0.10 → 10% chance of Type I Error

These values reflect the strictness of the test. As α increases, the likelihood of a Type I Error increases, making the test more prone to false positives.

2.4 – Type II Error

Type II Error (β) occurs when the null hypothesis is not rejected even though it is actually false. In other words, it is a false negative. The probability of committing a Type II Error is influenced by several factors: the sample size, the effect size (difference between true mean and hypothesized mean), the significance level (α), and the variability in the data.

The power of a test is defined as 1 − β and represents the probability of correctly rejecting a false null hypothesis. A higher power indicates a more sensitive test.

Step 1 – Definitions (for n = 30)

# Hypothesized population mean (H0)
mu0 <- pop_mean

# Assumed true population mean (alternative)
mu <- mu0 + 15000      # you can adjust this shift

# Standard deviation from full population
sigma <- sd(values)

# Effect size (difference between true mean and hypothesized mean)
delta <- abs(mu - mu0)

# Sample size
n <- 30

To compute Type II Error and power, we must assume a specific alternative population mean (μ) that differs from the null hypothesis mean (μ₀). A shift of 15,000 was chosen to represent a moderate, realistic difference in production values, enabling a meaningful power analysis.

Step 2-Power & β for α = 0.05 (n = 30)

# Significance level for first case
alpha <- 0.01

library(pwr)

# Compute power using z-test approximation
pwr_result_01 <- pwr.norm.test(d = delta / sigma, n = n, sig.level = alpha, alternative = "two.sided")
power_01 <- pwr_result_01$power
beta_01 <- 1 - power_01
power_01
## [1] 0.0249
beta_01
## [1] 0.975

Step 3-Power & β for α = 0.05 (n = 30)

alpha <- 0.05

pwr_result_05 <- pwr.norm.test(d = delta / sigma, n = n, sig.level = alpha, alternative = "two.sided")
power_05 <- pwr_result_05$power
beta_05 <- 1 - power_05
power_05
## [1] 0.0924
beta_05
## [1] 0.908

Step 4 – Power & β for α = 0.10 (n = 30)

alpha <- 0.10

pwr_result_10 <- pwr.norm.test(d = delta / sigma, n = n, sig.level = alpha, alternative = "two.sided")
power_10 <- pwr_result_10$power
beta_10 <- 1 - power_10
power_10
## [1] 0.161
beta_10
## [1] 0.839

Power and Type II Error for Different Alpha Levels

beta_df <- data.frame(
  `Alpha (α)` = c(0.01, 0.05, 0.10),
  `Power (1 - β)` = round(c(power_01, power_05, power_10), 4),
  `Type II Error (β)` = round(c(beta_01, beta_05, beta_10), 4)
)

knitr::kable(beta_df, caption = "Power and Type II Error for Different Alpha Levels")
Power and Type II Error for Different Alpha Levels
Alpha..α. Power..1…β. Type.II.Error..β.
0.01 0.025 0.975
0.05 0.092 0.908
0.10 0.161 0.839

The table above summarizes the calculated power and Type II Error values for three different significance levels. As expected:

  • When αlpha = 0.01, the test is more conservative, resulting in a lower power and higher Type II Error.
  • As αlpha increases, the power of the test also increases, since the acceptance region becomes narrower, making it easier to reject the null hypothesis when it is false.

This highlights the trade-off between Type I and Type II errors in hypothesis testing.

Step 6 – Recalculate for n = 40

n <- 40   #Increased sample size
alpha <- 0.01

pwr_result_01_n40 <- pwr.norm.test(d = delta / sigma, n = n, sig.level = alpha, alternative = "two.sided")
power_01_n40 <- pwr_result_01_n40$power
beta_01_n40 <- 1 - power_01_n40
alpha <- 0.05

pwr_result_05_n40 <- pwr.norm.test(d = delta / sigma, n = n, sig.level = alpha, alternative = "two.sided")
power_05_n40 <- pwr_result_05_n40$power
beta_05_n40 <- 1 - power_05_n40
alpha <- 0.10

pwr_result_10_n40 <- pwr.norm.test(d = delta / sigma, n = n, sig.level = alpha, alternative = "two.sided")
power_10_n40 <- pwr_result_10_n40$power
beta_10_n40 <- 1 - power_10_n40

Step 7 – Comparison Table

beta_compare <- data.frame(
  `Alpha(α)` = c(0.01, 0.05, 0.10),
  `Power(n=30)` = round(c(power_01, power_05, power_10), 4),
  `Power(n=40)` = round(c(power_01_n40, power_05_n40, power_10_n40), 4),
  `Beta(n=30)` = round(c(beta_01, beta_05, beta_10), 4),
  `Beta(n=40)` = round(c(beta_01_n40, beta_05_n40, beta_10_n40), 4)
)

knitr::kable(beta_compare, caption = "Comparison of Power and Type II Error for n = 30 and n = 40")
Comparison of Power and Type II Error for n = 30 and n = 40
Alpha.α. Power.n.30. Power.n.40. Beta.n.30. Beta.n.40.
0.01 0.025 0.030 0.975 0.970
0.05 0.092 0.107 0.908 0.893
0.10 0.161 0.181 0.839 0.819

The comparison table shows that increasing the sample size from 30 to 40 improves the power of the test across all significance levels. As the sample size increases, the standard error decreases, making the test more sensitive to differences and reducing the probability of Type II Error. This confirms that larger samples lead to more reliable and conclusive hypothesis tests.

Section 3 – Exploratory Data Analysis

In this section, an exploratory data analysis (EDA) was conducted to investigate the relationship between two major agricultural products in Argentina: Maize and Wheat. According to the project guidelines, data spanning the years 1961 to 2023 was used, focusing specifically on the variable Gross Production Value (constant 2014-2016 thousand I$).

Both data sets were obtained separately and then merged based on the Year variable to enable a side-by-side comparison of annual production values. This alignment ensures data consistency and allows for accurate correlation analysis between the two crop types.

3.1-Investigate the relationships between Wheat & Maize with Correlation Coefficent & Scatter Plot

Data Filtering & Reshaping

# Wheat data
wheat_data <- read.csv("~/Downloads/FAOSTAT_data_en_6-13-2025 (4).csv")

# Maize data
maize_data <- read.csv("~/Downloads/FAOSTAT_data_en_6-13-2025 (3).csv")
# Rename columns for clarity
wheat_data <- wheat_data[, c("Year", "Value")]
maize_data <- maize_data[, c("Year", "Value")]
colnames(wheat_data)[2] <- "Wheat"
colnames(maize_data)[2] <- "Maize"

# Merge by Year
combined_df <- merge(maize_data, wheat_data, by = "Year")

Correlation Coefficient Calculation

correlation <- cor(combined_df$Maize, combined_df$Wheat, use = "complete.obs")

# Printing the result
correlation
## [1] 0.706

The correlation coefficient between maize and wheat production in Argentina between 1961 and 2023 was calculated using the cor() function. The resulting value was approximately 0.706, indicating a moderately strong positive linear relationship between the two variables.

Scatterplot and Regression Line of Maize vs. Wheat

# Plot with regression line
plot(combined_df$Maize, combined_df$Wheat,
     main = "Scatterplot with Regression Line: Maize vs Wheat Production",
     xlab = "Maize Production Value (thousand I$)",
     ylab = "Wheat Production Value (thousand I$)",
     pch = 19, col = "orange")

# Add regression line
abline(lm(Wheat ~ Maize, data = combined_df), col = "purple", lwd = 2)

The scatter plot above visualizes the relationship between maize and wheat production in Argentina from 1961 to 2023. The upward pattern of points supports the positive correlation found in the numerical analysis, suggesting that these two crops tend to vary together over time.

A linear regression line was added to the scatter plot to highlight the direction and strength of the relationship between maize and wheat production. The positive slope confirms that as maize production increases, wheat production tends to increase as well.

3.2- Simple Linear Regression Model

model <- lm(Wheat ~ Maize, data = combined_df)
summary(model)
## 
## Call:
## lm(formula = Wheat ~ Maize, data = combined_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1442506  -526554   -31463   580943  1274266 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.85e+06   1.39e+05   13.30  < 2e-16 ***
## Maize       2.36e-01   3.04e-02    7.78  1.1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 713000 on 61 degrees of freedom
## Multiple R-squared:  0.498,  Adjusted R-squared:  0.49 
## F-statistic: 60.5 on 1 and 61 DF,  p-value: 1.06e-10

A simple linear regression model was fitted with wheat production as the dependent variable and maize production as the independent variable. The regression equation is of the form:

Wheat = β₀ + β₁ × Maize

According to the model summary: - Intercept (β₀): 1.846^{6} - Slope (β₁): 0.236

This suggests that for every additional unit increase in maize production value, wheat production increases by approximately 0.236 thousand I$ on average.

In conclusion, Simple linear equation is : Wheat=1,850,000+0.236×Maize

3.3- Residual Analysis

3.3.1 – Calculating Residual Analysis

residuals <- resid(model) ; residuals
##        1        2        3        4        5        6        7        8 
##  -720276  -743739    64366   566877  -650187  -700479  -516054  -797796 
##        9       10       11       12       13       14       15       16 
##  -508875 -1124747 -1028619  -253054  -752464  -901676  -181614   619194 
##       17       18       19       20       21       22       23       24 
##  -900871  -300260  -288752  -260965  -475609  1274266   832297   942921 
##       25       26       27       28       29       30       31       32 
##  -339813  -346051  -129295  -246359   302997   279291   403274    14444 
##       33       34       35       36       37       38       39       40 
##     4322   -31463   314536   -84955  1231761   808931   497780  1024050 
##       41       42       43       44       45       46       47       48 
##  1249626  1110082   377012   928803  1004969   481892   595009  1014490 
##       49       50       51       52       53       54       55       56 
##  -453291  -785309   830823   626361 -1442506 -1238794  -150442 -1053069 
##       57       58       59       60       61       62       63 
##   164651   478873    66667    69012  -537053   600686  -835826
fitted <- fitted(model) ;fitted
##       1       2       3       4       5       6       7       8       9      10 
## 2076132 2093674 2052901 2099837 2089881 2179961 2249655 2157204 2171427 2289954 
##      11      12      13      14      15      16      17      18      19      20 
## 2316978 2124017 2306074 2315556 2211253 2123780 2239699 2306074 2258663 2149619 
##      21      22      23      24      25      26      27      28      29      30 
## 2457788 2301333 2272886 2296592 2410377 2419859 2284739 2282369 2078502 2102208 
##      31      32      33      34      35      36      37      38      39      40 
## 2210532 2353508 2363014 2337365 2386864 2344869 2582801 2764092 2486429 2641772 
##      41      42      43      44      45      46      47      48      49      50 
## 2574390 2543700 2559462 2555019 2817283 2531063 2877627 2890028 2468284 2920663 
##      51      52      53      54      55      56      57      58      59      60 
## 2974556 2851137 3368984 3414875 3449513 3732797 4191876 3906769 4541995 4614776 
##      61      62      63 
## 4715760 4645183 3809440

3.3.2 – Residuals vs Fitted Plot

plot(fitted, residuals,
     main = "Residuals vs Fitted Values",
     xlab = "Fitted Values",
     ylab = "Residuals",
     pch = 19, col = "darkred")
abline(h = 0, col = "blue", lty = 2)

3.3.3 – Histogram of Residuals

hist(residuals,
     breaks = 10,
     main = "Histogram of Residuals",
     xlab = "Residuals",
     col = "lightblue", border = "black")

3.3.4 – QQ Plot

qqnorm(residuals)
qqline(residuals, col = "red")

The residuals vs. fitted plot shows a fairly even spread around the horizontal axis, suggesting that the linearity and homoscedasticity assumptions of the regression model are reasonably met. There is no clear pattern, such as curvature or funnel shape, which would indicate a violation of assumptions.

The histogram of residuals is approximately symmetric, and the QQ plot aligns reasonably well with the theoretical normal line. Together, these plots indicate that the residuals are roughly normally distributed.

Overall, the residual diagnostics support the validity of the linear regression model.

3.4- Confidence and Prediction Intervals for Two Values of the Regressor (Maize)

# Define new and arbitrary Maize values
new_data <- data.frame(Maize = c(1000000, 2000000))

# Confidence Intervals (for mean response)
conf_int <- predict(model, newdata = new_data, interval = "confidence")

# Prediction Intervals (for individual observations)
pred_int <- predict(model, newdata = new_data, interval = "prediction")
conf_int
##       fit     lwr     upr
## 1 2082379 1847745 2317013
## 2 2318568 2117513 2519624
pred_int
##       fit    lwr     upr
## 1 2082379 636570 3528189
## 2 2318568 877827 3759310

Confidence and prediction intervals were computed for two hypothetical maize production values: 1,000,000 and 2,000,000 (thousand I$). The confidence interval estimates the range in which the average wheat production is expected to fall, while the prediction interval accounts for individual variation and provides a wider range to capture possible new data points.

As expected, the prediction intervals are wider than the confidence intervals due to added uncertainty from individual-level noise.

Plotting CI ve PI

# Merging predictions as a data frame
conf_df <- as.data.frame(conf_int)
pred_df <- as.data.frame(pred_int)
new_data$Fitted <- conf_df$fit
new_data$CI_Lower <- conf_df$lwr
new_data$CI_Upper <- conf_df$upr
new_data$PI_Lower <- pred_df$lwr
new_data$PI_Upper <- pred_df$upr

# Graph
plot(new_data$Maize, new_data$Fitted,
     ylim = range(c(new_data$PI_Lower, new_data$PI_Upper)),
     xlab = "Maize Production Value (thousand I$)",
     ylab = "Predicted Wheat Production (thousand I$)",
     main = "Confidence and Prediction Intervals",
     pch = 19, col = "blue")

# CI
arrows(new_data$Maize, new_data$CI_Lower,
       new_data$Maize, new_data$CI_Upper,
       length = 0.05, angle = 90, code = 3, col = "darkgreen", lwd = 2)

# PI 
arrows(new_data$Maize, new_data$PI_Lower,
       new_data$Maize, new_data$PI_Upper,
       length = 0.05, angle = 90, code = 3, col = "orange", lwd = 2)

legend("topleft",
       legend = c("Fitted", "Confidence Interval", "Prediction Interval"),
       col = c("blue", "darkgreen", "orange"), pch = c(19, NA, NA),
       lty = c(NA, 1, 1), lwd = 2, bty = "n")

Section 4 – Trend-based forecasting

In this section, a trend-based forecasting approach is employed to estimate future production values for grapes in Argentina. The dataset contains annual gross production values (constant 2014–2016 I$) from 1961 onward.

A simple linear regression model is used to capture the long-term trend in grape production over time. By fitting a linear model of the form:

Production_Value = β₀ + β₁ × Year + ε

we can project the expected production values for future years. This approach assumes that the historical trend in grape production will continue into the near future.

The model will be used to forecast grape production values for the years 2025, 2030, and 2035.

4.1 – Linear Regression Model

# Read grapes dataset
grapes_raw <- read.csv("~/Downloads/FAOSTAT_data_en_6-14-2025.csv")

# Filter grapes for Argentina and correct element
grapes_data <- subset(grapes_raw, Area == "Argentina" & 
                      Item == "Grapes" &
                      Element == "Gross Production Value (constant 2014-2016 thousand I$)")

# Sort and select columns
grapes_data <- grapes_data[order(grapes_data$Year), c("Year", "Value")]

# Fit linear regression model
grapes_model <- lm(Value ~ Year, data = grapes_data)

# View model summary
summary(grapes_model)
## 
## Call:
## lm(formula = Value ~ Year, data = grapes_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -708658 -266616  -70777  268332  862071 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21452016    5413012    3.96  0.00020 ***
## Year           -9637       2717   -3.55  0.00076 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 392000 on 61 degrees of freedom
## Multiple R-squared:  0.171,  Adjusted R-squared:  0.157 
## F-statistic: 12.6 on 1 and 61 DF,  p-value: 0.000757

4.1.1 – Linear Regression Equation

coeffs <- coef(grapes_model)
cat("Regression Equation: ŷ =", round(coeffs[1], 2), "+", round(coeffs[2], 2), "* Year")
## Regression Equation: ŷ = 21452016 + -9637 * Year

4.2 – Forecasting for 2025, 2030, 2035

# Create new data frame for future years
future_years <- data.frame(Year = c(2025, 2030, 2035))

# Predict grape production values
predicted_values <- predict(grapes_model, newdata = future_years)

# Combine results
future_years$Predicted_Value <- round(predicted_values, 2)

# Display forecast
future_years
##   Year Predicted_Value
## 1 2025         1936267
## 2 2030         1888080
## 3 2035         1839893

4.2.1 – Manuel Calculation/Forecasting for 2025, 2030, 2035

To verify the predictions generated by the linear model, we can manually apply the regression equation to forecast grape production values for selected future years. The regression equation obtained from the model summary is as follows:

\[ \hat{y} = 21,\!452,\!016 - 9637 \times \text{Year} \]

Using this equation:

  • For the year 2025: \[ \hat{y}_{2025} = 21,\!452,\!016 - 9637 \times 2025 = 1,\!939,\!191 \]

  • For the year 2030: \[ \hat{y}_{2030} = 21,\!452,\!016 - 9637 \times 2030 = 1,\!891,\!006 \]

  • For the year 2035: \[ \hat{y}_{2035} = 21,\!452,\!016 - 9637 \times 2035 = 1,\!842,\!821 \]

These hand-calculated forecasts are consistent with the values produced by the predict() function, thereby validating the application of the trend-based model.

4.3 – Plotting Historical Data + Forecast Points

# Plot historical data
plot(grapes_data$Year, grapes_data$Value, type = "p",
     main = "Trend-Based Forecasting of Grape Production",
     xlab = "Year", ylab = "Production Value (thousand I$)",
     col = "darkgreen", pch = 16, xlim = c(1960, 2040))

# Add regression line (historical trend)
abline(grapes_model, col = "blue", lwd = 2)

# Add forecast points for 2025, 2030, 2035
points(future_years$Year, future_years$Predicted_Value, 
       col = "red", pch = 17, cex = 1.2)

# Add legend
legend("topleft",
       legend = c("Historical Data", "Forecast", "Trend Line"),
       col = c("darkgreen", "red", "blue"),
       pch = c(16, 17, NA), lty = c(NA, NA, 1), lwd = c(NA, NA, 2))

The plot above displays the historical grape production values from 1961 to 2023 in Argentina, along with a fitted linear trend line. Forecasted production values for the years 2025, 2030, and 2035 are also shown in red. These estimates are based on the continuation of the historical trend.

Conclusion

This project explored the agricultural production trends in Argentina using a dataset derived from FAOSTAT. We began with descriptive statistical analysis of apple production, constructing frequency tables and visualizations to summarize distribution characteristics.

Inferential statistics were then applied to perform hypothesis testing at various significance levels, enabling us to assess whether the sample mean significantly differed from the population mean. Confidence intervals were computed, and the statistical power of the tests was evaluated, including both Type I and Type II error analysis.

In the exploratory data analysis section, we investigated relationships between the production values of wheat and maize. The analysis revealed a moderately strong positive correlation, and a simple linear regression model was fitted to quantify this relationship. Residual diagnostics supported the model’s assumptions, and confidence/prediction intervals were estimated for selected input values.

Finally, a trend-based forecasting model was developed for grape production. A linear regression model was used to project future values for the years 2025, 2030, and 2035. Despite some variation in historical data, the trend model provided interpretable and consistent predictions.

Overall, this project applied key concepts from descriptive statistics, inferential analysis, and linear modeling to generate actionable insights from agricultural production data.

References

  1. R Studio Lab Manual – Marmara University IE2152 Statistics Course Materials
    Used for base R syntax, sampling methods, t-test procedures, and power analysis code structures in Part 1 and Part 2.

  2. Custom Code – Authored by the project owner
    Frequency distributions, descriptive summary tables, and data handling operations were implemented manually using custom R scripts.

  3. OpenAI ChatGPT . Technical assistance and code generation.
    AI-generated code was selectively used in Part 3 and Part 4 to assist with plotting residuals, visualizing regression lines, and formatting confidence/prediction intervals. All AI-supported outputs were actively reviewed, interpreted, and tested by the author to ensure compatibility with the dataset and accuracy of the statistical logic.
    No unverified or blindly copied output was used; model assumptions and graph structures were evaluated based on project context.

  4. FAOSTAT Dataset. Food and Agriculture Organization of the United Nations.
    https://www.fao.org/faostat/
    The dataset used in this study contains agricultural production indicators from 1961 to 2023 for Argentina.