function (Nsim = 1000) 
{
    treatment = as.factor(fertilizerdata$treatment)
    row = as.factor(fertilizerdata$row)
    response = fertilizerdata$response
    out.lm = lm(response ~ treatment + row)
    out = anova(out.lm)
    F.stat0 = out[4][[1]][1]
    F.stat = NULL
    for (i in 1:Nsim) {
        treatment.star = as.factor(c(sample(1:6), sample(1:6), 
            sample(1:6), sample(1:6)))
        out.lm = lm(response ~ treatment.star + row)
        out = anova(out.lm)
        F.stat = c(F.stat, out[4][[1]][1])
    }
    M = max(c(F.stat, 1.2 * F.stat0))
    hist(F.stat, nclass = 100, probability = T, xlim = c(0, M), 
        xlab = "randomization F-statistic", main = "Randomization Distribution (F-test for Treatment Effect)")
    abline(v = F.stat0, lwd = 2, col = "red")
    x = seq(0, M, 0.01)
    y = df(x, 5, 15)
    lines(x, y, col = "blue")
    text(1.02 * F.stat0, 0.2, paste("p-value =", format(signif(mean(F.stat >= 
        F.stat0), 5))), adj = 0, srt = 90)
    text(1.5, 0.7, expression("with" ~ F[5 ~ "," ~ 15] ~ "-density superimposed"), 
        adj = 0)
    text(1.5, 0.6, substitute(N[sim] == Nsim0, list(Nsim0 = Nsim)), 
        adj = 0)
    F.stat0
}