R Notebook

Install required packages

install.packages(c(“agricolae”, “car”, “ggplot2”, “dplyr”, “tidyr”))

Load libraries

library(gtsummary)
library(dplyr)
library(agricolae)
library(car)
library(ggplot2)
library(tidyr)

Load the built-in npk dataset

data(npk, package = "datasets")
head(npk)

str(npk)

'data.frame':   24 obs. of  5 variables:
 $ block: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ N    : Factor w/ 2 levels "0","1": 1 2 1 2 2 2 1 1 1 2 ...
 $ P    : Factor w/ 2 levels "0","1": 2 2 1 1 1 2 1 2 2 2 ...
 $ K    : Factor w/ 2 levels "0","1": 2 1 1 2 1 2 2 1 1 2 ...
 $ yield: num  49.5 62.8 46.8 57 59.8 58.5 55.5 56 62.8 55.8 ...

Check the structure

glimpse(npk)

Rows: 24
Columns: 5
$ block <fct> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6
$ N     <fct> 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0
$ P     <fct> 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0
$ K     <fct> 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0
$ yield <dbl> 49.5, 62.8, 46.8, 57.0, 59.8, 58.5, 55.5, 56.0, 62.8, 55.8, 69.5, 55.0, 62.0, 48.8, 45.5, 44.2, 52.0, 51.5, 49.8, 48.8, 5…

Data Exploration and Visualization

Basic summary statistics

summary(npk)

 block N      P      K          yield      
 1:4   0:12   0:12   0:12   Min.   :44.20  
 2:4   1:12   1:12   1:12   1st Qu.:49.73  
 3:4                        Median :55.65  
 4:4                        Mean   :54.88  
 5:4                        3rd Qu.:58.62  
 6:4                        Max.   :69.50  

Interpretation

It is fertilizer treatment data with yield summary statistics.

block, N, P, K → These are treatment factors:

block → Replication/block number (e.g., 1:4 means block 1 of 4).

N, P, K → Amounts of Nitrogen, Phosphorus, and Potassium applied.

yield → Crop yield obtained under each treatment.

The part you pasted under yield looks like the summary statistics generated by R (using summary(yield)), which gives:

Min. : 44.20 → Lowest yield recorded

1st Qu.: 49.73 → 25th percentile

Median : 55.65 → Middle value

Mean : 54.88 → Average yield

3rd Qu.: 58.62 → 75th percentile

Max. : 69.50 → Highest yield recorded

So in short:

👉 The crop yield ranged from 44.2 to 69.5, with an average around 55.

👉 Fertilizer treatments clearly influenced yield—some combinations produced significantly higher yields than others.

👉 Since the mean and median are close, yield distribution is fairly balanced, not skewed by a few extreme values.

Check for missing values

sum(is.na(npk))

[1] 0

Visualize yield distribution

ggplot(npk, aes(x = yield)) +
  geom_histogram(bins = 15, fill = "lightblue", color = "black") +
  ggtitle("Distribution of Yield")

Boxplot by treatment

ggplot(npk, aes(x = N, y = yield, fill = N)) +
  geom_boxplot() +
  ggtitle("Yield by Nitrogen Treatment")

Interaction plot

ggplot(npk, aes(x = N, y = yield, color = P, group = P)) +
  stat_summary(fun = mean, geom = "line", size = 1) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  ggtitle("Interaction between N and P on Yield")

Statistical Analysis

ANOVA Analysis

# Two-way ANOVA

anova_model <- aov(yield ~ N * P * K, data = npk)
summary(anova_model)

            Df Sum Sq Mean Sq F value Pr(>F)  
N            1  189.3  189.28   6.161 0.0245 *
P            1    8.4    8.40   0.273 0.6082  
K            1   95.2   95.20   3.099 0.0975 .
N:P          1   21.3   21.28   0.693 0.4175  
N:K          1   33.1   33.14   1.078 0.3145  
P:K          1    0.5    0.48   0.016 0.9019  
N:P:K        1   37.0   37.00   1.204 0.2887  
Residuals   16  491.6   30.72                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Check ANOVA assumptions

par(mfrow = c(2, 2))
plot(anova_model)
par(mfrow = c(1, 1))

Post-hoc tests

TukeyHSD(anova_model)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = yield ~ N * P * K, data = npk)

$N
        diff       lwr      upr     p adj
1-0 5.616667 0.8195759 10.41376 0.0245421

$P
         diff       lwr      upr     p adj
1-0 -1.183333 -5.980424 3.613757 0.6081875

$K
         diff       lwr       upr     p adj
1-0 -3.983333 -8.780424 0.8137574 0.0974577

$`N:P`
             diff        lwr       upr     p adj
1:0-0:0  7.500000  -1.655822 16.655822 0.1294203
0:1-0:0  0.700000  -8.455822  9.855822 0.9961506
1:1-0:0  4.433333  -4.722489 13.589156 0.5257140
0:1-1:0 -6.800000 -15.955822  2.355822 0.1873570
1:1-1:0 -3.066667 -12.222489  6.089156 0.7743737
1:1-0:1  3.733333  -5.422489 12.889156 0.6554324

$`N:K`
             diff        lwr        upr     p adj
1:0-0:0  7.966667  -1.189156 17.1224889 0.0999349
0:1-0:0 -1.633333 -10.789156  7.5224889 0.9554188
1:1-0:0  1.633333  -7.522489 10.7891555 0.9554188
0:1-1:0 -9.600000 -18.755822 -0.4441778 0.0382331
1:1-1:0 -6.333333 -15.489156  2.8224889 0.2364506
1:1-0:1  3.266667  -5.889156 12.4224889 0.7399770

$`P:K`
             diff       lwr      upr     p adj
1:0-0:0 -1.466667 -10.62249 7.689156 0.9670135
0:1-0:0 -4.266667 -13.42249 4.889156 0.5562955
1:1-0:0 -5.166667 -14.32249 3.989156 0.3986648
0:1-1:0 -2.800000 -11.95582 6.355822 0.8176375
1:1-1:0 -3.700000 -12.85582 5.455822 0.6615998
1:1-0:1 -0.900000 -10.05582 8.255822 0.9919305

$`N:P:K`
                    diff        lwr       upr     p adj
1:0:0-0:0:0  12.33333333  -3.335528 28.002195 0.1841773
0:1:0-0:0:0   2.90000000 -12.768862 18.568862 0.9975401
1:1:0-0:0:0   6.50000000  -9.168862 22.168862 0.8280891
0:0:1-0:0:0   0.56666667 -15.102195 16.235528 1.0000000
1:0:1-0:0:0   3.23333333 -12.435528 18.902195 0.9952141
0:1:1-0:0:0  -0.93333333 -16.602195 14.735528 0.9999987
1:1:1-0:0:0   2.93333333 -12.735528 18.602195 0.9973597
0:1:0-1:0:0  -9.43333333 -25.102195  6.235528 0.4627120
1:1:0-1:0:0  -5.83333333 -21.502195  9.835528 0.8902999
0:0:1-1:0:0 -11.76666667 -27.435528  3.902195 0.2249895
1:0:1-1:0:0  -9.10000000 -24.768862  6.568862 0.5045183
0:1:1-1:0:0 -13.26666667 -28.935528  2.402195 0.1303264
1:1:1-1:0:0  -9.40000000 -25.068862  6.268862 0.4668300
1:1:0-0:1:0   3.60000000 -12.068862 19.268862 0.9909550
0:0:1-0:1:0  -2.33333333 -18.002195 13.335528 0.9993808
1:0:1-0:1:0   0.33333333 -15.335528 16.002195 1.0000000
0:1:1-0:1:0  -3.83333333 -19.502195 11.835528 0.9870208
1:1:1-0:1:0   0.03333333 -15.635528 15.702195 1.0000000
0:0:1-1:1:0  -5.93333333 -21.602195  9.735528 0.8819049
1:0:1-1:1:0  -3.26666667 -18.935528 12.402195 0.9949094
0:1:1-1:1:0  -7.43333333 -23.102195  8.235528 0.7204567
1:1:1-1:1:0  -3.56666667 -19.235528 12.102195 0.9914323
1:0:1-0:0:1   2.66666667 -13.002195 18.335528 0.9985452
0:1:1-0:0:1  -1.50000000 -17.168862 14.168862 0.9999671
1:1:1-0:0:1   2.36666667 -13.302195 18.035528 0.9993214
0:1:1-1:0:1  -4.16666667 -19.835528 11.502195 0.9793323
1:1:1-1:0:1  -0.30000000 -15.968862 15.368862 1.0000000
1:1:1-0:1:1   3.86666667 -11.802195 19.535528 0.9863677

Linear Models

Multiple regression

lm_model <- lm(yield ~ N + P + K + N:P + N:K + P:K, data = npk)
summary(lm_model)

Call:
lm(formula = yield ~ N + P + K + N:P + N:K + P:K, data = npk)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8917 -3.2875  0.4583  3.4000  9.7083 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  52.6750     3.0114  17.492 2.64e-12 ***
N1            9.8500     3.9429   2.498    0.023 *  
P1            0.4167     3.9429   0.106    0.917    
K1           -1.9167     3.9429  -0.486    0.633    
N1:P1        -3.7667     4.5529  -0.827    0.420    
N1:K1        -4.7000     4.5529  -1.032    0.316    
P1:K1         0.5667     4.5529   0.124    0.902    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.576 on 17 degrees of freedom
Multiple R-squared:  0.3968,    Adjusted R-squared:  0.184 
F-statistic: 1.864 on 6 and 17 DF,  p-value: 0.146

Check model diagnostics

plot(lm_model)

Creating My Own NPK Dataset

Generate synthetic NPK data

set.seed(123)
n <- 100

synthetic_npk <- data.frame(
  block = factor(rep(1:4, each = n/4)),
  N = factor(sample(c("0", "1"), n, replace = TRUE)),
  P = factor(sample(c("0", "1"), n, replace = TRUE)),
  K = factor(sample(c("0", "1"), n, replace = TRUE)),
  yield = rnorm(n, mean = 50, sd = 10) +
    ifelse(sample(c("0", "1"), n, replace = TRUE) == "1", 5, 0) +
    ifelse(sample(c("0", "1"), n, replace = TRUE) == "1", 3, 0) +
    ifelse(sample(c("0", "1"), n, replace = TRUE) == "1", 2, 0)
)

head(synthetic_npk)

Advanced Analysis

Interaction Effects

Extract coefficients

coefficients <- coef(interaction_model)
print(coefficients)

(Intercept)          N1          P1          K1       N1:P1       N1:K1       P1:K1    N1:P1:K1 
 51.4333333  12.3333333   2.9000000   0.5666667  -8.7333333  -9.6666667  -4.4000000   9.9333333 

Calculate effect sizes

library(effects)
effect_plot <- allEffects(interaction_model)
plot(effect_plot)

Treatment Combinations Analysis

Create treatment combinations

npk <- npk %>%
  mutate(treatment = interaction(N, P, K, sep = "_"))

Compare all treatment combinations

treatment_aov <- aov(yield ~ treatment, data = npk)
summary(treatment_aov)

            Df Sum Sq Mean Sq F value Pr(>F)
treatment    7  384.8   54.97   1.789  0.159
Residuals   16  491.6   30.72               

Post-hoc comparison

tukey <- TukeyHSD(treatment_aov)
print(tukey)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = yield ~ treatment, data = npk)

$treatment
                    diff        lwr       upr     p adj
1_0_0-0_0_0  12.33333333  -3.335528 28.002195 0.1841773
0_1_0-0_0_0   2.90000000 -12.768862 18.568862 0.9975401
1_1_0-0_0_0   6.50000000  -9.168862 22.168862 0.8280891
0_0_1-0_0_0   0.56666667 -15.102195 16.235528 1.0000000
1_0_1-0_0_0   3.23333333 -12.435528 18.902195 0.9952141
0_1_1-0_0_0  -0.93333333 -16.602195 14.735528 0.9999987
1_1_1-0_0_0   2.93333333 -12.735528 18.602195 0.9973597
0_1_0-1_0_0  -9.43333333 -25.102195  6.235528 0.4627120
1_1_0-1_0_0  -5.83333333 -21.502195  9.835528 0.8902999
0_0_1-1_0_0 -11.76666667 -27.435528  3.902195 0.2249895
1_0_1-1_0_0  -9.10000000 -24.768862  6.568862 0.5045183
0_1_1-1_0_0 -13.26666667 -28.935528  2.402195 0.1303264
1_1_1-1_0_0  -9.40000000 -25.068862  6.268862 0.4668300
1_1_0-0_1_0   3.60000000 -12.068862 19.268862 0.9909550
0_0_1-0_1_0  -2.33333333 -18.002195 13.335528 0.9993808
1_0_1-0_1_0   0.33333333 -15.335528 16.002195 1.0000000
0_1_1-0_1_0  -3.83333333 -19.502195 11.835528 0.9870208
1_1_1-0_1_0   0.03333333 -15.635528 15.702195 1.0000000
0_0_1-1_1_0  -5.93333333 -21.602195  9.735528 0.8819049
1_0_1-1_1_0  -3.26666667 -18.935528 12.402195 0.9949094
0_1_1-1_1_0  -7.43333333 -23.102195  8.235528 0.7204567
1_1_1-1_1_0  -3.56666667 -19.235528 12.102195 0.9914323
1_0_1-0_0_1   2.66666667 -13.002195 18.335528 0.9985452
0_1_1-0_0_1  -1.50000000 -17.168862 14.168862 0.9999671
1_1_1-0_0_1   2.36666667 -13.302195 18.035528 0.9993214
0_1_1-1_0_1  -4.16666667 -19.835528 11.502195 0.9793323
1_1_1-1_0_1  -0.30000000 -15.968862 15.368862 1.0000000
1_1_1-0_1_1   3.86666667 -11.802195 19.535528 0.9863677

Data Visualization Practice

Create comprehensive visualizations

1. Faceted plots

ggplot(npk, aes(x = block, y = yield, fill = N)) +
  geom_boxplot() +
  facet_grid(P ~ K) +
  ggtitle("Yield by Block, N, P, and K Treatments")

2. Interaction heatmap

interaction_means <- npk %>%
  group_by(N, P, K) %>%
  summarise(mean_yield = mean(yield), .groups = "drop")

ggplot(interaction_means, aes(x = N, y = P, fill = mean_yield)) +
  geom_tile() +
  facet_wrap(~ K) +
  scale_fill_gradient(low = "blue", high = "red") +
  ggtitle("Mean Yield by Treatment Combinations")

Practice Exercises

Exercise 1: Which treatment combination gives maximum yield?

max_yield_treatment <- npk %>%
  group_by(N, P, K) %>%
  summarise(mean_yield = mean(yield), .groups = "drop") %>%
  arrange(desc(mean_yield)) %>%
  head(1)

print(max_yield_treatment)

NA

Exercise 2: Test if there’s a significant difference between blocks

block_anova <- aov(yield ~ block, data = npk)
summary(block_anova)

            Df Sum Sq Mean Sq F value Pr(>F)  
block        5  343.3   68.66   2.318 0.0861 .
Residuals   18  533.1   29.62                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Exercise 3: Create a predictive model
set.seed(123)
train_indices <- sample(1:nrow(npk), 0.7 * nrow(npk))
train_data <- npk[train_indices, ]
test_data <- npk[-train_indices, ]

pred_model <- lm(yield ~ N + P + K + N:P, data = train_data)
predictions <- predict(pred_model, newdata = test_data)

Calculate RMSE

rmse <- sqrt(mean((test_data$yield - predictions)^2))
print(paste("RMSE:", round(rmse, 2)))

[1] "RMSE: 8.4"

Exporting Results

Generate report

sink("npk_analysis_report.txt")
cat("NPK Dataset Analysis Report\n")
cat("===========================\n\n")
cat("Dataset Summary:\n")
print(summary(npk))
cat("\nANOVA Results:\n")
print(summary(anova_model))
sink()