knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

# Install once if needed:
# install.packages(c("agricolae","gt","broom","corrplot","car"))

library(agricolae)  # sweetpotato and corn datasets
library(gt)          # nice tables
library(dplyr)        # select() used with tidy() output
library(broom)          # tidy() turns test/model output into a data frame
library(corrplot)         # correlogram
library(car)                # vif() for multicollinearity check

1 Dataset

Module 7 uses corn, a different built-in agricolae dataset — a completely randomized design comparing 4 corn-growing methods. It has 34 plots and, importantly, two numeric variables recorded on the same plots:

method — growing method (1 to 4)
observation — yield per plot
rx — a second numeric measurement taken on the same plot

Having two numeric variables on the same units is what lets us do real (not simulated) correlation, regression, and paired testing in Module 7.

data(sweetpotato)
data(corn)

str(sweetpotato)

## 'data.frame':    12 obs. of  2 variables:
##  $ virus: Factor w/ 4 levels "cc","fc","ff",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ yield: num  28.5 21.7 23 14.9 10.6 13.1 41.8 39.2 28 38.2 ...

str(corn)

## 'data.frame':    34 obs. of  3 variables:
##  $ method     : int  1 1 1 1 1 1 1 1 1 2 ...
##  $ observation: int  83 91 94 89 89 96 91 92 90 91 ...
##  $ rx         : num  11 23 28.5 17 17 31.5 23 26 19.5 23 ...

2 MODULE 7: Inferential Statistics

Descriptive stats describe the sample. Inferential stats let us draw conclusions about the wider population, and say how confident we are in those conclusions. Everything below uses the corn dataset.

corn$method <- factor(corn$method)

2.1 Confidence Interval (CI)

Definition. A range of plausible values for a population parameter (here, the true mean observation), at a chosen confidence level (usually 95%).

What it’s for. Shows how precise our estimate of the mean really is.

Note. Data is normal, or sample size large enough; observations independent.

ci_result <- t.test(corn$observation)

tidy(ci_result) |>
  select(estimate, conf.low, conf.high) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 2) |>
  tab_header(title = "95% Confidence Interval for Mean Observation (Yield)")

estimate	conf.low	conf.high
95% Confidence Interval for Mean Observation (Yield)
87.88	85.64	90.13

Interpretation. We are 95% confident the true mean yield falls within conf.low and conf.high.

2.2 Hypothesis Testing

2.2.1 One-Sample t-test

Definition. Tests if a sample mean is significantly different from a fixed value. Used to compare the mean value of a sample with a constant value.

What it’s for. Example: is the average observation different from a benchmark of 90?

Assumptions. Randomly selected. Normally distributed or n > 30, independent observations, no significant outliers.

t_one <- t.test(corn$observation, mu = 90)

tidy(t_one) |>
  select(estimate, statistic, p.value, conf.low, conf.high) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 3) |>
  tab_header(title = "One-Sample t-test: Observation vs. Benchmark (90)")

estimate	statistic	p.value	conf.low	conf.high
One-Sample t-test: Observation vs. Benchmark (90)
87.882	−1.921	0.063	85.640	90.125

Interpretation. At 5% level of significance, we can say there is no significant difference from the claim value of 90.

2.2.2 Two Independent Sample t-test

Definition. Compares the means of two independent groups. Used to compare the mean values of two independent samples, to determine whether they are drawn from populations with equal means.

What it’s for. Example: comparing observation between method 1 and method 2.

Assumptions. Normality of data in each group; similar variances; independent observations.

# Keep only method 1 and 2
two_groups <- subset(corn, method %in% c(1, 2))
two_groups$method <- factor(two_groups$method)

# Check equal variance assumption first
var_check <- var.test(observation ~ method, data = two_groups)

tidy(var_check) |>
  select(estimate, statistic, p.value, conf.low, conf.high) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 3) |>
  tab_header(title = "Equal Variance Check (F-test): Method 1 vs. Method 2")

estimate	statistic	p.value	conf.low	conf.high
Equal Variance Check (F-test): Method 1 vs. Method 2
0.931	0.931	0.930	0.227	4.055

# Two-sample t-test
t_two <- t.test(observation ~ method, data = two_groups)

tidy(t_two) |>
  select(estimate1, estimate2, statistic, p.value, conf.low, conf.high) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 3) |>
  tab_header(title = "Two-Sample t-test: Observation, Method 1 vs. Method 2")

estimate1	estimate2	statistic	p.value	conf.low	conf.high
Two-Sample t-test: Observation, Method 1 vs. Method 2
90.556	86.400	2.439	0.026	0.560	7.751

Interpretation. For the F-test: p.value ≥ 0.05 means the equal variance assumption is reasonable.

For the t-test: p.value < 0.05 means the two methods have significantly different mean yields;

At 5% level of significance, we can conclude that there is a significant difference of the mean yields of Method 1 and Method 2.

2.2.3 Paired t-test

Definition. Used to compare the mean values for two samples, where each value in one sample corresponds to a particular value in the other sample.

The following are the scores of eight students before and after a review. At 5% level of significance, test score is different before and after the review.

Student	1	2	3	4	5	6	7	8
*Before*	77	74	82	73	87	68	66	80
*After*	72	68	76	68	84	68	61	76

Assumptions. The differences between pairs are roughly normal.

Before <- c(77, 74, 82, 73, 87, 68, 66, 80)
After <- c(72, 68, 76, 68, 84, 68, 61, 76)


t_paired <- t.test(Before, After, paired = TRUE) 

tidy(t_paired) |>    
  select(estimate, statistic, p.value, conf.low, conf.high) |>   
  gt() |>    
  fmt_number(columns = everything(), decimals = 3) |>  
  tab_header(title = "Paired t-test: Before vs. After")

estimate	statistic	p.value	conf.low	conf.high
Paired t-test: Before vs. After
4.250	6.065	0.001	2.593	5.907

Interpretation. At 5% level of significance, we can say that is a significant difference between Before and after of the exam.

2.2.4 Chi-Square Test of Independence

Definition. The chi-square test of association (sometimes called the chi-square test of independence) helps to determine whether two or more categorical variables are associated.

Ho: The row variable and column variable are not related

Ha: The row variable and column variable are related

What it’s for. Example: is “high vs. low yield” related to growing method?

Assumptions. Sample data are randomly selected. For every cell in the contingency table, the expected frequency is at least 5.

corn$yield_level <- ifelse(corn$observation >= median(corn$observation),
                            "High", "Low")

chi_table <- table(corn$method, corn$yield_level)
chi_table

##    
##     High Low
##   1    8   1
##   2    4   6
##   3    7   0
##   4    0   8

chi_result <- chisq.test(chi_table)

tidy(chi_result) |>
  gt() |>
  fmt_number(columns = c(statistic, p.value), decimals = 4) |>
  tab_header(title = "Chi-Square Test: Method vs. Yield Level")

statistic	p.value	parameter	method
Chi-Square Test: Method vs. Yield Level
20.6598	0.0001	3	Pearson's Chi-squared test

chisq.test(chi_table)$expected

##    
##         High      Low
##   1 5.029412 3.970588
##   2 5.588235 4.411765
##   3 3.911765 3.088235
##   4 4.470588 3.529412

Interpretation. p.value < 0.05 suggests yield level and growing method are related (not independent).

2.2.5 Nonparametric Tests

Used when data isn’t normally distributed, or the sample is small.

Wilcoxon test (alternative to two-sample t-test):

wilcox_result <- wilcox.test(observation ~ method, data = two_groups)

tidy(wilcox_result) |>
  gt() |>
  fmt_number(columns = c(statistic, p.value), decimals = 4) |>
  tab_header(title = "Wilcoxon Rank-Sum Test: Method 1 vs. Method 2")

statistic	p.value	method	alternative
Wilcoxon Rank-Sum Test: Method 1 vs. Method 2
70.5000	0.0395	Wilcoxon rank sum test with continuity correction	two.sided

2.3 Correlation

2.3.1 Pearson Correlation with Significance Test

Definition.

A hypothesis test of correlation determines whether a correlation is statistically significant.
The null hypothesis for the test is that the population correlation is equal to zero, meaning that there is no correlation between the variables.
The alternative hypothesis is that the population correlation is not equal to zero, meaning that there is some correlation between the variables.
You can also perform a one-sided test, where the alternative hypothesis is either that the population correlation is greater than zero (the variables are positively correlated) or that the population correlation is less than zero (the variables are negatively correlated). • You can perform a test of the correlation between two variables with the cor.test function:

cor.test(dataset$var1, dataset$var2)
By default, R performs a test of the Pearson’s correlation. If you would prefer to test the Spearman’s correlation, set the method argument to “spearman”:

cor.test(dataset$var1, dataset$var2, method="spearman")
By default, R performs a two-sided test, but you can adjust this by setting the alternative argument to “less” or “greater” as required:

cor.test(dataset$var1, dataset$var2, alternative="greater")
The output includes a 95% confidence interval for the correlation estimate. To adjust the size of this interval, use the conf.level argument:

cor.test(dataset$var1, dataset$var2, conf.level=0.99)

Value	Interval	Interpretation
Very weak	0.00 - 0.19	Almost no linear relationship
Weak	0.20 – 0.39	Small but noticeable relationship
Moderate	0.40 – 0.59	Clear, moderate relationship
Strong	0.60 – 0.79	Strong, consistent relationship
Very strong	0.80 – 1.00	Very strong, almost perfect relationship

What it’s for. Does rx relate to observation (yield)?

Assumptions. Both variables numeric; relationship is linear; no extreme outliers.

# ============================================================
# CORRELATION (using iris instead of corn)
# ============================================================

data(iris)
str(iris)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

## Pearson Correlation with Significance Test ---------------

# Definition: measures the strength/direction of a LINEAR relationship
# between two numeric variables (-1 to +1).
# What it's for: does Sepal.Length relate to Petal.Length?
# Assumptions: both numeric; linear relationship; no extreme outliers.

pearson_result <- cor.test(iris$Sepal.Length, iris$Petal.Length, method = "pearson")

tidy(pearson_result) |>
  select(estimate, statistic, p.value, conf.low, conf.high) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 3) |>
  tab_header(title = "Pearson Correlation: Sepal.Length vs. Petal.Length")

estimate	statistic	p.value	conf.low	conf.high
Pearson Correlation: Sepal.Length vs. Petal.Length
0.872	21.646	0.000	0.827	0.906

Interpretation. estimate (r) is about 0.87 — a strong positive relationship: as Sepal.Length increases, Petal.Length tends to increase too. p.value is essentially 0 (< 0.001), so this relationship is highly statistically significant.

## Correlogram -------------------------------------------------

# Definition: a picture of the correlation between several variables at once.
# What it's for: quickly spot which variables are strongly related,
# before building a regression model.

iris_num <- iris[, c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")]
corr_matrix <- cor(iris_num)

corrplot(corr_matrix, method = "color", addCoef.col = "white")

Interpretation. You’ll see Petal.Length and Petal.Width correlate very strongly (~0.96) with each other, Sepal.Length correlates strongly with both petal measurements (~0.87 and ~0.82), while Sepal.Width correlates weakly, even negatively, with the others. This screening step is exactly why we picked Sepal.Length as a predictor of Petal.Length for regression below — and why Sepal.Width alone wouldn’t be a good simple-regression predictor.

## Spearman Rank Correlation ------------------------------------

spearman_result <- cor.test(iris$Sepal.Length, iris$Petal.Length, method = "spearman")

tidy(spearman_result) |>
  select(estimate, statistic, p.value) |>
  gt() |>
  fmt_number(columns = c(estimate, p.value), decimals = 3) |>
  tab_header(title = "Spearman Rank Correlation: Sepal.Length vs. Petal.Length")

estimate	statistic	p.value
Spearman Rank Correlation: Sepal.Length vs. Petal.Length
0.882	66429.35	0.000

Interpretation. Spearman’s rho comes out close to the Pearson r (~0.88), suggesting the relationship is both linear and monotonic — the two methods agree.

2.4 Regression

2.4.1 Simple Linear Regression

Definition. Used to predict the value of a dependent variable based on the value of at least one independent variable.

Used to explain the impact of changes in an independent variable on the dependent variable.

Components of Regression Analysis:

Dependent/ Response variable – the variable we wish to explain
Independent/ Predictor / Explanatory variable – the variable used to explain the dependent variable

To build a simple linear regression model with an explanatory variable named var1 and a response variable named resp, use the command:

lm(resp~var1, dataset)
You don’t need to specify an intercept term (or constant term) in your model because R includes one automatically. To build a model without an intercept term, use the command:

lm(resp~-1+var1, dataset)
You can save all the output to an object name using this command:

modelname<-lm(resp~var1, dataset)

Assumptions:

Linearity: The relationship between X and the mean of Y is linear.
Homoscedasticity: The variance of residual is the same for any value of X.
Independence: Observations are independent of each other.
Normality: For any fixed value of X, Y is normally distributed.

# ============================================================
# SIMPLE LINEAR REGRESSION
# ============================================================

# Definition: models one numeric outcome (Y) using one numeric predictor (X).
# What it's for: predicting Petal.Length from Sepal.Length.
# Assumptions (LINE): Linearity, Independence, Normality of residuals,
# Equal variance of residuals.

model_slr <- lm(Petal.Length ~ Sepal.Length, data = iris)

tidy(model_slr) |>
  gt() |>
  fmt_number(columns = c(estimate, std.error, statistic, p.value), decimals = 4) |>
  tab_header(title = "Simple Linear Regression: Petal.Length ~ Sepal.Length")

term	estimate	std.error	statistic	p.value
Simple Linear Regression: Petal.Length ~ Sepal.Length
(Intercept)	−7.1014	0.5067	−14.0161	0.0000
Sepal.Length	1.8584	0.0859	21.6460	0.0000

glance(model_slr) |>
  select(r.squared, adj.r.squared, sigma, statistic, p.value) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 4) |>
  tab_header(title = "Model Fit: Petal.Length ~ Sepal.Length")

r.squared	adj.r.squared	sigma	statistic	p.value
Model Fit: Petal.Length ~ Sepal.Length
0.7600	0.7583	0.8678	468.5502	0.0000

Expected output (standard values for this well-known pairing):

Intercept ≈ -7.10
Sepal.Length coefficient (slope) ≈ 1.86
R-squared ≈ 0.76

Interpretation — Regression Equation:

Each 1 cm increase in Sepal.Length is associated with about a 1.86 cm increase in predicted Petal.Length. The intercept (-7.10) has no real biological meaning on its own here (a flower with 0 cm sepal length doesn’t exist) — it’s just the line’s starting point mathematically.

Interpretation — Test of Significance of the Slope:

The p.value for the Sepal.Length row tests H₀: slope = 0 (no relationship) vs. H₁: slope ≠ 0. With p.value < 0.001, we reject H₀ — Sepal.Length is a statistically significant predictor of Petal.Length.

Interpretation — Coefficient of Determination (R²):

R² ≈ 0.76 means about 76% of the variation in Petal.Length is explained by Sepal.Length alone. The remaining ~24% is due to other factors not in the model (e.g., species differences).

## Diagnostic plots to evaluate predictive ability -----------

par(mfrow = c(2, 2))
plot(model_slr)

par(mfrow = c(1, 1))

# Simple scatter with fitted line
plot(iris$Sepal.Length, iris$Petal.Length,
     xlab = "Sepal Length", ylab = "Petal Length",
     main = "Simple Linear Regression Fit")
abline(model_slr, col = "red", lwd = 2)

Interpretation — Plots. The scatter plot with fitted line shows a clear upward trend, though points fan out somewhat, especially at higher Sepal.Length values — a hint that Species might explain some of that spread (worth flagging as a segue into the multiple regression below). In the diagnostic plots: Residuals vs. Fitted should show a rough horizontal band (if there’s a curve, linearity is questionable); Q-Q plot points should hug the diagonal (normality of residuals); Scale-Location should show a fairly flat spread (equal variance); and no point should stand out with extreme Cook’s distance (undue influence).

2.4.2 Multiple Linear Regression

# ============================================================
# MULTIPLE LINEAR REGRESSION
# ============================================================

# Definition: same idea, but with more than one predictor.
# What it's for: predicting Petal.Length from Sepal.Length,
# Sepal.Width, AND Petal.Width together.
# Assumptions: same as simple regression, plus low multicollinearity
# among predictors (check with VIF).

model_mlr <- lm(Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width, data = iris)

tidy(model_mlr) |>
  gt() |>
  fmt_number(columns = c(estimate, std.error, statistic, p.value), decimals = 4) |>
  tab_header(title = "Multiple Linear Regression: Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width")

term	estimate	std.error	statistic	p.value
Multiple Linear Regression: Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width
(Intercept)	−0.2627	0.2974	−0.8833	0.3785
Sepal.Length	0.7291	0.0583	12.5025	0.0000
Sepal.Width	−0.6460	0.0685	−9.4312	0.0000
Petal.Width	1.4468	0.0676	21.3987	0.0000

glance(model_mlr) |>
  select(r.squared, adj.r.squared, sigma, statistic, p.value) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 4) |>
  tab_header(title = "Model Fit: Multiple Regression")

r.squared	adj.r.squared	sigma	statistic	p.value
Model Fit: Multiple Regression
0.9680	0.9674	0.3190	1,472.7262	0.0000

# Multicollinearity check
vif(model_mlr)

## Sepal.Length  Sepal.Width  Petal.Width 
##     3.415733     1.305515     3.889961

Interpretation — Regression Equation (fill in your actual coefficients after running):

Each b coefficient is the expected change in Petal.Length per 1-unit increase in that predictor, holding the other two predictors constant — this is the key difference from simple regression’s interpretation.

Interpretation — Test of Significance of Each Slope: Check each predictor’s p.value individually. A predictor with p.value < 0.05 significantly contributes to explaining Petal.Length even after accounting for the other two predictors; one with a high p-value may be redundant once the others are in the model (watch especially for Sepal.Length, since it’s correlated with Petal.Width, which can weaken its individual significance in the presence of the other predictor — this is exactly what vif() checks for).

Interpretation — R² and Comparison to Simple Regression: Compare adj.r.squared here to the simple model’s R² (0.76). Adding Sepal.Width and Petal.Width typically pushes R² noticeably higher (multiple regression using all three usually explains 95%+ of the variance in Petal.Length in this dataset) — a strong illustration for trainees of why adding relevant predictors improves prediction, as long as vif() values stay reasonably low (below ~5) confirming the predictors aren’t just duplicating each other’s information.

Statistical Analysis in R

Module 7: Inferential Statistics

Rice Data Analytics