Disclaimer: I had look through my notes from School to remind me of a few things, I also used Claude AI for some assistance. So Nick, Matt and Sanjeev. If you see something wrong. Please let me know.

Design: This was the first field Gilbron and I collected data from, before we started digging in pairs, so it is set up differently from 230 and 263. It is one variety, comparing whole vs cut seed, and the digs were taken at random spots in the field instead of in pairs. With no pairs there are no blocks, so I analyzed it as two independent groups (6 whole digs, 6 cut digs), the way a completely randomized design (CRD) is analyzed, using a Welch two-sample t-test.

âš  Important - please read: why this field is analyzed differently

  • The digs were taken at random spots, not in pairs, so there are no blocks and I could not use a paired t-test like in 230 and 263.
  • Instead, the 6 whole digs are compared with the 6 cut digs as two groups (Welch t-test).
  • The GPS shows the whole and cut digs ended up in different parts of the field (nearest whole and cut digs about 440 m apart). So any difference - or lack of one - could partly be due to location, not just seed type.

1 Make sure the data is clean

1.1 Load the data

xl <- "Field548_data.xlsx"            # keep this file in the same folder as this .Rmd
if (!file.exists(xl)) stop("Can't find ", xl, " in ", getwd(), " - put it in the same folder as this .Rmd")
d <- as.data.frame(readxl::read_excel(xl, sheet = "Data"))
d$Treatment <- factor(d$Treatment, levels = c("Whole", "Cut"))
str(d)
## 'data.frame':    12 obs. of  18 variables:
##  $ Dig          : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ Treatment    : Factor w/ 2 levels "Whole","Cut": 1 1 1 2 2 1 1 2 2 1 ...
##  $ Stems        : num  19 18 28 22 17 18 19 20 23 19 ...
##  $ n_10oz       : num  16 25 16 13 23 24 13 23 17 24 ...
##  $ lbs_10oz     : num  13.2 20.2 13.2 9.8 18 18 9.8 17.8 13 20.4 ...
##  $ n_2s         : num  4 4 1 3 2 3 3 2 1 1 ...
##  $ lbs_2s       : num  1.8 2.8 1.4 1.96 0.12 ...
##  $ n_4.5to10oz  : num  44 32 33 49 37 39 34 42 20 37 ...
##  $ lbs_4.5to10oz: num  19.6 15 15.2 21.6 15.4 17.2 15.4 18.4 18.6 16 ...
##  $ n_undersize  : num  32 5 17 17 12 19 10 14 14 20 ...
##  $ lbs_undersize: num  6.8 1 3.4 4 2.6 3.4 2 4 3.4 4 ...
##  $ n_process    : num  1 3 12 2 3 1 3 1 1 1 ...
##  $ lbs_process  : num  0.8 1.2 6.4 1.5 1.47 ...
##  $ Rev_Jumbo    : num  7648 11703 7648 5678 10429 ...
##  $ Rev_2s       : num  505 785.6 392.8 549.9 33.7 ...
##  $ Rev_ASize    : num  7483 5726 5803 8246 5879 ...
##  $ Rev_UndProc  : num  759 220 979 549 407 ...
##  $ Revenue      : num  16395 18435 14822 15023 16748 ...

This is just to make sure the data loads correctly and looks clean.

1.2 Let me check for any missing values

colSums(is.na(d))
##           Dig     Treatment         Stems        n_10oz      lbs_10oz 
##             0             0             0             0             0 
##          n_2s        lbs_2s   n_4.5to10oz lbs_4.5to10oz   n_undersize 
##             0             0             0             0             0 
## lbs_undersize     n_process   lbs_process     Rev_Jumbo        Rev_2s 
##             0             0             0             0             0 
##     Rev_ASize   Rev_UndProc       Revenue 
##             0             0             0

No missing values here, data looks good.

1.3 Design check: 6 digs of each treatment

table(d$Treatment)
## 
## Whole   Cut 
##     6     6
stopifnot(all(table(d$Treatment) == 6))   # stops knitting if a treatment label is mistyped

Each treatment must show 6. stopifnot() is a safety check: if “Whole” or “Cut” were mistyped in the Excel file, knitting would stop here. And that is why I use this function.

1.4 Weight check: average ounces per tuber

oz_per <- function(lbs, n) ifelse(n > 0, 16 * lbs / n, NA)   # lb per tuber -> oz; NA if no tubers
chk <- data.frame(
  Dig    = d$Dig,
  Jumbo  = oz_per(d$lbs_10oz,        d$n_10oz),          # expect >= 10 oz
  A_size = oz_per(d$`lbs_4.5to10oz`, d$`n_4.5to10oz`),   # expect 4.5-10 oz
  Under  = oz_per(d$lbs_undersize,   d$n_undersize),     # expect < 4.5 oz (culls can be larger)
  No2s   = oz_per(d$lbs_2s,          d$n_2s))
knitr::kable(round(chk, 1), row.names = FALSE, caption = "Average oz per tuber by size class (NA = none in that class)")
Average oz per tuber by size class (NA = none in that class)
Dig Jumbo A_size Under No2s
1 13.2 7.1 3.4 7.2
2 12.9 7.5 3.2 11.2
3 13.2 7.4 3.2 22.4
4 12.1 7.1 3.8 10.5
5 12.5 6.7 3.5 1.0
6 12.0 7.1 2.9 6.4
7 12.1 7.2 3.2 5.3
8 12.4 7.0 4.6 16.8
9 12.2 14.9 3.9 7.0
10 13.6 6.9 3.2 2.1
11 13.2 7.0 3.2 5.4
12 11.8 7.1 3.2 7.6 
zero_wt <- with(d, (n_2s > 0 & lbs_2s == 0) | (n_process > 0 & lbs_process == 0))
if (any(zero_wt)) cat("WARNING - tubers counted but 0 lb, digs:", d$Dig[zero_wt], "\n")

Each average should fall inside its size class. This check found the original data errors (#2s and process weights entered in ounces instead of pounds, and two weights left at 0), which were fixed from the field notes. One value still doesn’t fit: dig 9’s A size averages 14.9 oz (18.6 lb, 20 tubers). That matches the notes, so the tuber count was probably mis-tallied in the field; the weight is fine, so yield and revenue are not affected.

1.5 Calculate yield and % jumbos

row_ft <- 10; space_in <- 34                              # one 10-ft row at 34-in spacing
cwt_per_lb <- 43560 / (row_ft * space_in / 12) / 100      # lbs per dig -> cwt/acre

d$total_lbs   <- with(d, lbs_10oz + lbs_2s + `lbs_4.5to10oz` + lbs_undersize + lbs_process)
d$Yield_cwt   <- d$total_lbs * cwt_per_lb
d$Pct_jumbo   <- 100 * d$lbs_10oz / d$total_lbs
# Revenue columns come from the Excel file (Data sheet; prices on the Prices sheet)
stopifnot(all(abs(d$Rev_Jumbo + d$Rev_2s + d$Rev_ASize + d$Rev_UndProc - d$Revenue) < 1))

# Main results were chosen before looking at the data; supporting traits explain WHY
primary   <- c("Revenue", "Yield_cwt")
secondary <- c("Rev_Jumbo", "Rev_ASize", "Rev_2s", "Rev_UndProc", "Pct_jumbo", "Stems")
traits    <- c(primary, secondary)

gW <- d[d$Treatment == "Whole", ]
gC <- d[d$Treatment == "Cut", ]

1.6 Means by treatment

means <- aggregate(d[traits], by = list(Treatment = d$Treatment), FUN = mean)
means[traits] <- round(means[traits], 1)
knitr::kable(means, row.names = FALSE, caption = "Means by treatment (revenue in $/acre, yield in cwt/acre, stems per 10 ft)")
Means by treatment (revenue in $/acre, yield in cwt/acre, stems per 10 ft)
Treatment Revenue Yield_cwt Rev_Jumbo Rev_ASize Rev_2s Rev_UndProc Pct_jumbo Stems
Whole 16327.0 596.8 9154.1 6261.0 389.6 522.3 40.4 20.2
Cut 16782.4 616.3 8362.3 7711.7 253.9 454.5 36.1 19.0

These are just averages. They show how big the gaps are, but not whether they are real or could be field variation - that is what the t-tests in section 3 answer.

2 Check assumptions

With two independent groups, the analysis compares the 6 whole digs with the 6 cut digs, so the assumptions are:

  • each group is roughly normal - judged mainly from the Q-Q plot; Shapiro-Wilk is shown but has little power with 6 digs per group;
  • no single dig dominates - outlier check;
  • the residuals vs predicted plot shows no funnel or curve (course Transform handout).

Equal variance between the two groups is not needed for the Welch test, but I show it (F-test) because the ANOVA in section 4 does assume it.

assump <- data.frame()
for (v in traits) {
  fit <- aov(as.formula(paste(v, "~ Treatment")), data = d)
  sw  <- shapiro.test(residuals(fit))              # normality of the 12 residuals
  ft  <- var.test(gW[[v]], gC[[v]])                # F-test: equal variances?
  out <- c(gW$Dig[gW[[v]] %in% boxplot.stats(gW[[v]])$out],   # digs beyond 1.5 x IQR in their group
           gC$Dig[gC[[v]] %in% boxplot.stats(gC[[v]])$out])
  assump <- rbind(assump, data.frame(
    Trait        = v,
    Role         = ifelse(v %in% primary, "Main", "Supporting"),
    Shapiro_P    = round(sw$p.value, 3),
    SD_Whole     = round(sd(gW[[v]]), 1),
    SD_Cut       = round(sd(gC[[v]]), 1),
    EqualVar_P   = round(ft$p.value, 3),
    Outlier_digs = ifelse(length(out) == 0, "none", paste(out, collapse = ","))))
}
knitr::kable(assump, row.names = FALSE, caption = "Assumption checks (P > 0.05 = no evidence against normality / equal variance)")
Assumption checks (P > 0.05 = no evidence against normality / equal variance)
Trait Role Shapiro_P SD_Whole SD_Cut EqualVar_P Outlier_digs
Revenue Main 0.163 2448.9 1548.7 0.337 none
Yield_cwt Main 0.197 69.3 53.9 0.594 7
Rev_Jumbo Supporting 0.400 2523.8 1794.2 0.472 none
Rev_ASize Supporting 0.557 670.9 1391.4 0.135 none
Rev_2s Supporting 0.320 248.7 247.2 0.990 none
Rev_UndProc Supporting 0.493 286.5 73.9 0.010 none
Pct_jumbo Supporting 0.321 8.9 7.9 0.807 none
Stems Supporting 0.397 3.9 4.9 0.617 3 
fig1 <- function() {
  base_par(); par(mfrow = c(2, 2), mar = c(4.2, 4.6, 3, 1), oma = c(1.5, 0, 0, 0))
  for (v in primary) {
    fit <- aov(as.formula(paste(v, "~ Treatment")), data = d)
    r   <- rstandard(fit)                                         # standardized residuals
    lab <- paste("dig", d$Dig)
    out <- abs(r) > 2                                             # flag |residual| > 2
    ttl <- ifelse(v == "Revenue", "Total revenue", "Total yield")
    q <- qqnorm(r, plot.it = FALSE)
    plot(q$x, q$y, xlim = range(q$x) * 1.25, ylim = range(q$y) * 1.15, pch = 21, cex = 1.4,
         bg = col_v[as.character(d$Treatment)], col = "white",
         xlab = "Theoretical quantiles", ylab = "Standardized residuals",
         main = paste("Q-Q plot -", ttl), font.main = 1, col.main = ink,
         panel.first = grid(col = grid_c))
    qqline(r, col = ink, lty = 2)
    if (any(out)) text(q$x[out], q$y[out], lab[out], pos = 4, cex = 0.8, col = ink)
    xr <- range(fitted(fit)); xr[2] <- xr[2] + diff(xr) * 0.12
    plot(fitted(fit), r, xlim = xr, ylim = range(c(r, -2.2, 2.2)), pch = 21, cex = 1.4,
         bg = col_v[as.character(d$Treatment)], col = "white",
         xlab = "Predicted", ylab = "Standardized residuals",
         main = paste("Residuals vs predicted -", ttl), font.main = 1, col.main = ink,
         panel.first = grid(col = grid_c))
    abline(h = 0, lty = 2, col = ink); abline(h = c(-2, 2), lty = 3, col = "grey60")
    if (any(out)) text(fitted(fit)[out], r[out], lab[out], pos = 4, cex = 0.8, col = ink)
    legend("topleft", names(col_v), pt.bg = col_v, pch = 21, col = "white", pt.cex = 1.3,
           bty = "n", cex = 0.85, text.col = ink)
  }
  mtext("Points outside +/-2 (dotted lines) are possible outliers (labelled by dig number)",
        side = 1, outer = TRUE, line = 0.3, cex = 0.75, col = "grey40")
}
show_and_save(fig1, "Fig1_QQ_residuals.png", 9.6, 8.4)
Figure 1. Q-Q plots and residuals vs predicted for the main traits.
Figure 1. Q-Q plots and residuals vs predicted for the main traits.

Take-home message for Figure 1:

  • Q-Q plot: the points follow the dashed line, so the data are roughly normal (Shapiro P = 0.163 for revenue and 0.197 for yield).
  • Residuals vs predicted: the points are scattered around 0 with no funnel or curve, so the model fits. (They sit in two columns because there are only two treatments.)
  • Dig 7 (whole) is at the edge: it had one of the lowest jumbo weights (9.8 lb). It is real data, and removing it would only move the whole-seed average closer to cut, so it is kept.
  • Dig 3 has the most stems (28) and shows up as an outlier for stems in the table. It matches the field notes, so it is kept.

3 Pick the right test

The design picks the test: two independent groups (no pairs) -> Welch two-sample t-test, chosen before seeing the data. It compares the two group averages without assuming the groups vary equally (two-sample handout; JMP: Fit Y by X -> t Test, unequal variances). No LSD is needed - with only 2 treatments, the test itself is the mean comparison. The Wilcoxon rank-sum test is shown only as an extra check.

results <- data.frame()
for (v in traits) {
  tt <- t.test(gW[[v]], gC[[v]], var.equal = FALSE)                        # main test (Welch)
  wt <- suppressWarnings(wilcox.test(gW[[v]], gC[[v]], exact = FALSE))      # extra check only
  results <- rbind(results, data.frame(
    Trait       = v,
    Role        = ifelse(v %in% primary, "Main", "Supporting"),
    Whole       = round(mean(gW[[v]]), 1),
    Cut         = round(mean(gC[[v]]), 1),
    Difference  = round(mean(gW[[v]]) - mean(gC[[v]]), 1),
    CI95_low    = round(tt$conf.int[1], 1),
    CI95_high   = round(tt$conf.int[2], 1),
    t           = round(unname(tt$statistic), 2),
    P_Welch     = signif(tt$p.value, 3),
    Significant = ifelse(tt$p.value < alpha, "Yes", "No"),
    P_Wilcoxon  = signif(wt$p.value, 3)))
}
rownames(results) <- NULL
knitr::kable(results[, names(results) != "Role"], row.names = FALSE,
             caption = "Welch t-tests, Whole - Cut, 6 digs each (first two rows = main results)")
Welch t-tests, Whole - Cut, 6 digs each (first two rows = main results)
Trait Whole Cut Difference CI95_low CI95_high t P_Welch Significant P_Wilcoxon
Revenue 16327.0 16782.4 -455.4 -3158.2 2247.4 -0.38 0.7100 No 0.8100
Yield_cwt 596.8 616.3 -19.5 -100.0 61.0 -0.55 0.5980 No 0.9360
Rev_Jumbo 9154.1 8362.3 791.8 -2066.7 3650.4 0.63 0.5470 No 0.5720
Rev_ASize 6261.0 7711.7 -1450.7 -2933.3 31.9 -2.30 0.0539 No 0.0542
Rev_2s 389.6 253.9 135.7 -183.3 454.7 0.95 0.3660 No 0.4710
Rev_UndProc 522.3 454.5 67.9 -232.1 367.8 0.56 0.5960 No 0.8100
Pct_jumbo 40.4 36.1 4.3 -6.6 15.1 0.88 0.4000 No 0.4710
Stems 20.2 19.0 1.2 -4.6 6.9 0.46 0.6570 No 0.8080

3.1 What do the difference and confidence interval mean?

  • Difference = the whole-seed average minus the cut-seed average. It is our best single estimate of how much more (or less) whole seed gives per acre.
  • 95% confidence interval (CI) = the range of values the true difference could reasonably be, given only 6 digs per treatment. A narrow interval means we know the difference well; a wide one means we don’t. (Formally: if the trial were repeated many times, 95% of intervals built this way would contain the true difference.)
  • CI and P-value agree: if the interval does not include 0, then P < 0.05 and the difference is significant. The CI is more useful because it also shows how big the difference could be.

Examples from this field:

Trait Difference 95% CI What it means
Total revenue -$455/acre -$3,158 to +$2,247 Anywhere from cut earning about $3,200 more to whole earning about $2,200 more. The range includes 0, so it is not significant.
Total yield -20 cwt/acre -100 to +61 Anywhere from cut yielding 100 cwt more to whole yielding 61 cwt more. The range includes 0, so it is not significant - we can’t tell, which is not the same as “they yield the same”.

Take-home message: the difference tells you the size of the effect; the CI tells you how sure you can be about it.

4 ANOVA tables

Model for each trait: Trait ~ Treatment (CRD, one factor). There is no Pair line because this field has no pairs. With 2 treatments, the ANOVA P-value equals the regular (pooled) t-test, which assumes equal variance - so it can differ slightly from the Welch P-value.

4.1 Total revenue (main)

fit_rev <- aov(Revenue ~ Treatment, data = d)
summary(fit_rev)
##             Df   Sum Sq Mean Sq F value Pr(>F)
## Treatment    1   622121  622121   0.148  0.708
## Residuals   10 41977865 4197786
cv_rev <- 100 * sqrt(deviance(fit_rev) / df.residual(fit_rev)) / mean(d$Revenue)
cat("CV % =", round(cv_rev, 1))
## CV % = 12.4

Take-home message for the above ANOVA table:

  • Seed type did not affect revenue: Treatment P = 0.708 (Welch P = 0.710).
  • CV = 12.4%: a measure of how much the digs varied (lower = more consistent). Revenue varied a bit more than yield (CV = 10.2%).

4.2 Total yield (main)

fit_yld <- aov(Yield_cwt ~ Treatment, data = d)
summary(fit_yld)
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1   1145    1145   0.297  0.598
## Residuals   10  38508    3851
cv_yld <- 100 * sqrt(deviance(fit_yld) / df.residual(fit_yld)) / mean(d$Yield_cwt)
cat("CV % =", round(cv_yld, 1))
## CV % = 10.2

Take-home message for the above ANOVA table:

  • Seed type did not affect yield: Treatment P = 0.598. Whole and cut were only 20 cwt/acre apart.
  • CV = 10.2%: the lowest of the three fields - the digs were fairly consistent.

4.3 Supporting traits

for (v in secondary) {
  cat("\n------------------------------------------------\nANOVA -", v, "\n")
  fit <- aov(as.formula(paste(v, "~ Treatment")), data = d)
  print(summary(fit))
  cat("CV % =", round(100 * sqrt(deviance(fit) / df.residual(fit)) / mean(d[[v]]), 1), "\n")
}
## 
## ------------------------------------------------
## ANOVA - Rev_Jumbo 
##             Df   Sum Sq Mean Sq F value Pr(>F)
## Treatment    1  1880883 1880883   0.392  0.545
## Residuals   10 47945163 4794516               
## CV % = 25 
## 
## ------------------------------------------------
## ANOVA - Rev_ASize 
##             Df   Sum Sq Mean Sq F value Pr(>F)  
## Treatment    1  6313697 6313697   5.292 0.0442 *
## Residuals   10 11930730 1193073                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## CV % = 15.6 
## 
## ------------------------------------------------
## ANOVA - Rev_2s 
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1  55218   55218   0.898  0.366
## Residuals   10 614946   61495               
## CV % = 77.1 
## 
## ------------------------------------------------
## ANOVA - Rev_UndProc 
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1  13812   13812   0.316  0.587
## Residuals   10 437766   43777               
## CV % = 42.8 
## 
## ------------------------------------------------
## ANOVA - Pct_jumbo 
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1   54.8   54.83   0.772    0.4
## Residuals   10  710.0   71.00               
## CV % = 22 
## 
## ------------------------------------------------
## ANOVA - Stems 
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1   4.08   4.083    0.21  0.657
## Residuals   10 194.83  19.483               
## CV % = 22.5

Take-home message for the supporting traits:

  • Nothing is significant with the Welch test.
  • A size revenue is the closest: cut was higher, Welch P = 0.054. The ANOVA above shows P = 0.044, but that assumes both groups vary equally, and cut varied about twice as much as whole here. So I report the Welch result: a borderline trend, not a significant difference.

5 Graphs

5.1 Individual digs: total revenue and yield

fig2 <- function() {
  base_par(); par(mfrow = c(1, 2), mar = c(3.5, 6, 3.5, 1.5), mgp = c(4, 0.6, 0), oma = c(1.5, 0, 0, 0))
  for (v in primary) {
    ttl <- ifelse(v == "Revenue", "Total revenue ($/acre)", "Total yield (cwt/acre)")
    p   <- results$P_Welch[results$Trait == v]
    yl  <- range(d[[v]]) * c(0.9, 1.08)
    plot(NA, xlim = c(0.4, 2.6), ylim = yl, xaxt = "n", yaxt = "n", xlab = "", ylab = ttl,
         main = paste0(ttl, "\n", pfmt(p), " ", stars(p)), font.main = 1, col.main = ink, cex.main = 0.95,
         panel.first = abline(h = pretty(yl), col = grid_c))
    axis(2, at = pretty(yl), labels = trimws(format(pretty(yl), big.mark = ",")), lwd = 0)
    axis(1, at = 1:2, labels = names(col_v), tick = FALSE)
    for (k in 1:2) {
      g  <- if (k == 1) gW else gC
      xs <- k + seq(-0.15, 0.15, length.out = nrow(g))[rank(g[[v]], ties.method = "first")]
      segments(k - 0.25, mean(g[[v]]), k + 0.25, mean(g[[v]]), lwd = 3, col = ink)
      points(xs, g[[v]], pch = 21, cex = 1.6, bg = col_v[k], col = "white", lwd = 1.5)
      text(xs, g[[v]], g$Dig, pos = ifelse(k == 1, 2, 4), cex = 0.65, col = "grey40")
    }
  }
  mtext("Dots = individual digs (dig number); black bars = treatment means",
        side = 1, outer = TRUE, line = 0.3, cex = 0.75, col = "grey40")
}
show_and_save(fig2, "Fig2_digs_yield_revenue.png", 9.6, 5.2)
Figure 2. Each dot is one dig; black bars are treatment means.
Figure 2. Each dot is one dig; black bars are treatment means.

Take-home message: - Revenue and yield: the whole and cut dots overlap almost completely, which is why neither difference is significant. - There are no lines here (unlike 230 and 263) because the digs were not taken in pairs.

5.2 Yield and revenue by size class

cls  <- c("Jumbo\n10oz+", "A size\n4.5-10oz", "#2s", "Undersize\n+ process")
ycol <- list(d$lbs_10oz, d$`lbs_4.5to10oz`, d$lbs_2s, d$lbs_undersize + d$lbs_process)
for (k in 1:4) d[[paste0("Yld_c", k)]] <- ycol[[k]] * cwt_per_lb   # yield per size class (cwt/acre)
yv <- paste0("Yld_c", 1:4); rv <- c("Rev_Jumbo", "Rev_ASize", "Rev_2s", "Rev_UndProc")
gW <- d[d$Treatment == "Whole", ]; gC <- d[d$Treatment == "Cut", ]   # refresh with new columns

size_panel <- function(vars, ylab, ttl) {
  m  <- rbind(sapply(vars, function(v) mean(gW[[v]])), sapply(vars, function(v) mean(gC[[v]])))
  se <- rbind(sapply(vars, function(v) sd(gW[[v]]) / sqrt(6)), sapply(vars, function(v) sd(gC[[v]]) / sqrt(6)))
  pv <- sapply(vars, function(v) t.test(gW[[v]], gC[[v]], var.equal = FALSE)$p.value)
  top <- max(m + se) * 1.35
  bp <- barplot(m, beside = TRUE, col = col_v, border = "white", space = c(0.08, 0.8),
                names.arg = cls, ylim = c(0, top), ylab = ylab, main = ttl, font.main = 1,
                col.main = ink, cex.names = 0.8, axes = FALSE,
                panel.first = abline(h = pretty(c(0, top)), col = grid_c))
  axis(2, at = pretty(c(0, top)), labels = trimws(format(pretty(c(0, top)), big.mark = ",")), lwd = 0)
  arrows(bp, m - se, bp, m + se, angle = 90, code = 3, length = 0.03, col = ink)
  text(colMeans(bp), apply(m + se, 2, max) + top * 0.05, stars(pv), col = ink, cex = 0.9)
  legend("topright", names(col_v), fill = col_v, border = "white", bty = "n", cex = 0.85, text.col = ink)
}
fig3 <- function() {
  base_par(); par(mfrow = c(1, 2), mar = c(4.5, 6, 3, 1), mgp = c(4, 0.6, 0), oma = c(1.5, 0, 0, 0))
  size_panel(yv, "Yield (cwt/acre)", "Yield by size class")
  size_panel(rv, "Revenue ($/acre)", "Revenue by size class")
  mtext("Bars = mean +/- SE (6 digs each). Welch t-test: * P < 0.05, ** P < 0.01, *** P < 0.001, ns = not significant",
        side = 1, outer = TRUE, line = 0.3, cex = 0.75, col = "grey40")
}
show_and_save(fig3, "Fig3_size_class_yield_revenue.png", 9.6, 5.4)
Figure 3. Mean +/- SE by size class; stars from Welch t-tests.
Figure 3. Mean +/- SE by size class; stars from Welch t-tests.

Take-home message: - Yield: not significantly different in any size class. - Size: whole seed leaned slightly toward jumbos, cut seed toward A size. - Revenue: the two leanings cancel out, so total revenue ends up about the same.

5.3 Revenue difference with 95% confidence intervals

fig4 <- function() {
  base_par(); par(mar = c(5, 9, 3, 6.5))
  lab4 <- c("Total revenue", "Jumbo 10oz+", "A size 4.5-10oz", "#2s", "Undersize + process")
  v4   <- c("Revenue", rv)
  ci   <- t(sapply(v4, function(v) { t <- t.test(gW[[v]], gC[[v]])
                                      c(unname(t$estimate[1] - t$estimate[2]), t$conf.int, t$p.value) }))
  yy <- rev(seq_along(v4))
  xl <- range(c(0, ci[, 2:3])) * 1.05
  plot(NA, xlim = xl, ylim = c(0.5, length(v4) + 0.5), yaxt = "n", xaxt = "n",
       xlab = "Whole minus Cut ($/acre)", ylab = "",
       main = "Revenue difference with 95% confidence interval", font.main = 1, col.main = ink,
       panel.first = abline(v = pretty(xl), col = grid_c))
  axis(1, at = pretty(xl), labels = trimws(format(pretty(xl), big.mark = ",")), lwd = 0)
  axis(2, at = yy[1], labels = lab4[1], lwd = 0, font = 2)
  axis(2, at = yy[-1], labels = lab4[-1], lwd = 0)
  abline(v = 0, col = ink, lwd = 1.2)
  sig <- ci[, 4] < alpha
  segments(ci[, 2], yy, ci[, 3], yy, lwd = 2, col = ifelse(sig, ink, "grey65"))
  points(ci[, 1], yy, pch = 21, cex = 1.6, bg = ifelse(ci[, 1] > 0, col_v[1], col_v[2]),
         col = "white", lwd = 1.5)
  for (i in seq_along(yy))
    mtext(paste(pfmt(ci[i, 4]), stars(ci[i, 4])), side = 4, at = yy[i], las = 1, line = 0.5,
          cex = 0.75, col = ifelse(sig[i], ink, "grey50"))
  mtext("Interval not crossing 0 = significant. Blue dot = Whole higher, orange = Cut higher",
        side = 1, line = 3.6, cex = 0.7, col = "grey40")
}
show_and_save(fig4, "Fig4_revenue_difference_CI.png", 8.8, 4.8)
Figure 4. Whole minus Cut revenue, with 95% confidence intervals.
Figure 4. Whole minus Cut revenue, with 95% confidence intervals.

Take-home message: - Every interval crosses 0, so no revenue category is significantly different. - A size comes closest (cut higher).

6 Summary

The story: whole and cut seed gave about the same yield and revenue in Field 548. Cut was 20 cwt/acre and $455/acre higher on average, but neither difference is significant (P = 0.598 and 0.710). The only hint of a difference is in tuber size: cut leaned toward A size (P = 0.054) and whole toward jumbos, and those cancel out in total revenue.

Methods sentence: Data were analyzed as a completely randomized design with two independent groups (Welch two-sample t-test), alpha = 0.05, in R. Normality and homogeneity of residuals were checked with Q-Q and residuals vs predicted plots, and equal variance with an F-test.

Limitations:

  • The whole and cut digs ended up in different parts of the field (nearest whole and cut digs about 440 m apart), so a seed effect can’t be fully separated from a location effect.
  • One field, one year, 6 digs per treatment - results are preliminary. Only a fairly large difference (roughly 20%) could have been detected.
  • Each dig is a single 10-ft row scaled to per acre.
  • Revenue depends on the $/cwt prices on the Prices sheet of Field548_data.xlsx.
write.csv(assump,  "Field548_assumptions.csv", row.names = FALSE)
write.csv(results, "Field548_results.csv",     row.names = FALSE)

Saved: Field548_results.csv, Field548_assumptions.csv, and Figures 1-4 as PNG files in this folder.