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 if the review is effective.

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

Kruskal-Wallis test (alternative to one-way ANOVA, compares more than two groups):

kw_result <- kruskal.test(observation ~ method, data = corn)

tidy(kw_result) |>
  gt() |>
  fmt_number(columns = c(statistic, p.value), decimals = 4) |>
  tab_header(title = "Kruskal-Wallis Test: Observation across All 4 Methods")

statistic	p.value	parameter	method
Kruskal-Wallis Test: Observation across All 4 Methods
25.6288	0.0000	3	Kruskal-Wallis rank sum test

Interpretation. Same rule for both: p.value < 0.05 means the groups’ medians are significantly different.

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.

pearson_result <- cor.test(corn$rx, corn$observation, 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: rx vs. Observation")

estimate	statistic	p.value	conf.low	conf.high
Pearson Correlation: rx vs. Observation
0.987	34.599	0.000	0.974	0.993

Interpretation. estimate is r. A positive value means as rx increases, observation tends to increase too. p.value < 0.05 means the relationship is statistically significant.

2.3.2 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.

# method is included here as a numeric code (1-4) purely to show a
# correlogram with more than 2 variables; treat it as illustrative,
# since method is really a categorical treatment, not a continuous scale
corn_num <- data.frame(
  observation = corn$observation,
  rx = corn$rx,
  method_code = as.numeric(as.character(corn$method))
)

corr_matrix <- cor(corn_num)

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

Interpretation. Darker boxes and coefficients closer to +1 or -1 mean a stronger relationship. observation and rx are the pair worth trusting here, since method_code isn’t truly continuous.

2.3.3 Spearman Rank Correlation

Definition. A nonparametric correlation, based on ranks instead of raw values. Captures relationships that go consistently up or down, even if not perfectly straight-line.

What it’s for. Used when data isn’t normal, or is ordinal.

Assumptions. Data at least ordinal; relationship is monotonic.

spearman_result <- cor.test(corn$rx, corn$observation, 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: rx vs. Observation")

estimate	statistic	p.value
Spearman Rank Correlation: rx vs. Observation
1.000	0	0.000

Interpretation. Same idea as Pearson, but based on ranks. If it differs a lot from the Pearson result, the relationship may not be a straight line.

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.

model1 <- lm(observation ~ rx, data = corn)

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

term	estimate	std.error	statistic	p.value
Simple Linear Regression: Observation ~ rx
(Intercept)	76.6983	0.3703	207.1248	0.0000
rx	0.6391	0.0185	34.5985	0.0000

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

r.squared	adj.r.squared	sigma	statistic	p.value
Model Fit: Observation ~ rx
0.9740	0.9732	1.0533	1,197.0583	0.0000

Interpretation. The rx coefficient tells you how much observation changes for each 1-unit increase in rx. r.squared tells you what percent of the variation in observation is explained by rx alone.

2.4.2 Multiple Linear Regression

Definition. Same idea, but with more than one predictor.

What it’s for. Predicting observation from both rx and method together.

Assumptions. Same as simple regression, plus: predictors shouldn’t be too strongly correlated with each other (multicollinearity).

model2 <- lm(observation ~ rx + method, data = corn)

tidy(model2) |>
  gt() |>
  fmt_number(columns = c(estimate, std.error, statistic, p.value), decimals = 4) |>
  tab_header(title = "Multiple Linear Regression: Observation ~ rx + Method")

term	estimate	std.error	statistic	p.value
Multiple Linear Regression: Observation ~ rx + Method
(Intercept)	76.5291	0.9648	79.3214	0.0000
rx	0.6424	0.0409	15.7207	0.0000
method2	0.0417	0.5721	0.0728	0.9424
method3	0.1875	0.6388	0.2936	0.7712
method4	0.2542	0.8776	0.2897	0.7741

glance(model2) |>
  select(r.squared, adj.r.squared, sigma, statistic, p.value) |>
  gt() |>
  fmt_number(columns = everything(), decimals = 4) |>
  tab_header(title = "Model Fit: Observation ~ rx + Method")

r.squared	adj.r.squared	sigma	statistic	p.value
Model Fit: Observation ~ rx + Method
0.9742	0.9706	1.1013	273.8022	0.0000

Interpretation. Each coefficient shows the effect of that predictor while holding the other predictors constant. Compare adj.r.squared here to the simple model’s r.squared — a higher value means adding method genuinely improved the model.

Statistical Analysis in R

Module 7: Inferential Statistics

Rice Data Analytics