function (Nsim = 1) 
{
    life = Battery$life
    temperature = as.factor(Battery$temp)
    MaterialType = as.factor(Battery$type)
    out.lm.full0 = lm(life ~ temperature * MaterialType)
    out0 = anova(out.lm.full0)
    F.temp0 = out0[[4]][1]
    F.type0 = out0[[4]][2]
    F.inter0 = out0[[4]][3]
    combs = rep(1, 36)
    F.temp = NULL
    F.type = NULL
    F.inter = NULL
    for (i in 1:Nsim) {
        life.star = sample(life)
        out.lm.full = lm(life.star ~ temperature * MaterialType)
        out = anova(out.lm.full)
        F.temp = c(F.temp, out[[4]][1])
        F.type = c(F.type, out[[4]][2])
        F.inter = c(F.inter, out[[4]][3])
    }
    M = max(c(F.temp, 1.06 * F.temp0, 35))
    hist(F.temp, xlim = c(0, M), ylim = c(0, 1), nclass = 50, 
        probability = T, xlab = "F for temperature", main = substitute("Randomization Reference Distribution," ~ 
            N[sim] == Nsim0, list(Nsim0 = Nsim)), col = c("orange", 
            "green"))
    abline(v = F.temp0, lwd = 2, col = "red")
    xx = seq(0, M, 0.01)
    yy = df(xx, 2, 27)
    lines(xx, yy, col = "blue", lwd = 2)
    text(0.98 * F.temp0, 0.6, substitute(F[temperature] == Ftemp, 
        list(Ftemp = signif(F.temp0, 4))), adj = 0.5, col = "red", 
        srt = -90)
    text(1.02 * F.temp0, 0.6, paste("F-distribution p-value =", 
        format(signif(out0[[5]][1], 4))), adj = 0.5, srt = 90, 
        col = "red")
    text(1.06 * F.temp0, 0.6, paste("Randomization p-value =", 
        format(signif(mean(F.temp >= F.temp0), 4))), adj = 0.5, 
        srt = 90, col = "red")
    readline("hit return\n")
    M = max(c(F.type, 1.1 * F.type0, 9))
    hist(F.type, xlim = c(0, M), ylim = c(0, 1), nclass = 50, 
        probability = T, xlab = "F for type", main = substitute("Randomization Reference Distribution," ~ 
            N[sim] == Nsim0, list(Nsim0 = Nsim)), col = c("orange", 
            "green"))
    abline(v = F.type0, lwd = 2, col = "red")
    xx = seq(0, M, 0.01)
    yy = df(xx, 2, 27)
    lines(xx, yy, col = "blue", lwd = 2)
    text(0.96 * F.type0, 0.6, substitute(F[type] == Ftype, list(Ftype = signif(F.type0, 
        4))), adj = 0.5, col = "red", srt = -90)
    text(1.03 * F.type0, 0.6, paste("F-distribution p-value =", 
        format(signif(out0[[5]][2], 4))), adj = 0.5, srt = 90, 
        col = "red")
    text(1.08 * F.type0, 0.6, paste("Randomization p-value =", 
        format(signif(mean(F.type >= F.type0), 4))), adj = 0.5, 
        srt = 90, col = "red")
    readline("hit return\n")
    M = max(c(F.inter, 1.08 * F.inter0, 8))
    M = 8
    hist(F.inter, xlim = c(0, M), ylim = c(0, 1), nclass = 50, 
        probability = T, xlab = "F for interactions", main = substitute("Randomization Reference Distribution," ~ 
            N[sim] == Nsim0, list(Nsim0 = Nsim)), col = c("orange", 
            "green"))
    abline(v = F.inter0, lwd = 2, col = "red")
    xx = seq(0, M, 0.01)
    yy = df(xx, 4, 27)
    lines(xx, yy, col = "blue", lwd = 2)
    text(0.96 * F.inter0, 0.6, substitute(F[interaction] == Finter, 
        list(Finter = signif(F.inter0, 4))), adj = 0.5, col = "red", 
        srt = -90)
    text(1.04 * F.inter0, 0.6, paste("F-distribution p-value =", 
        format(signif(out0[[5]][3], 4))), adj = 0.5, srt = 90, 
        col = "red")
    text(1.09 * F.inter0, 0.6, paste("Randomization p-value =", 
        format(signif(mean(F.inter >= F.inter0), 4))), adj = 0.5, 
        srt = 90, col = "red")
}