"psconf" <- function(type,numfac,kvec,numrep,idfac,ybarj,xbarj,xbar,
s2unbl,s2bl,s2anocov,b,conflev) {
# This function is based on research that is reported in
# Verrill and Kretschmann (2016), "A Reminder about Potentially
# Serious Problems with a Type of Blocked ANOVA Analysis," USDA Forest
# Products Laboratory draft research paper.
# This function calculates a confidence interval on a treatment mean
# from ANOVA/ANOCOV results obtained
# after a predictor sort sample allocation. Currently,
# as input, the function requires
# type = 0 if an analysis of covariance has not been performed
# = 1 if an analysis of covariance has been performed
# numfac is the number of factors (1 for a 1-way anova with reps,
# 2 for a two-way anova with reps, and so on).
# kvec is the vector of numbers of levels corresponding
# to the factors.
# numrep is the number of replicates.
# Thus, for example, for a 1-way anova that compares 5 treatments
# and has 10 replicates per treatment, numfac = 1,
# kvec = c(5), numrep = 10. For a 3 x 2 x 2 anova with 4 replicates
# per cell, numfac = 3, kvec = c(3,2,2), numrep = 4.
# idfac is the ID of the factor under consideration. The function will be
# calculating a confidence interval on the mean of one of its levels.
# ybarj is the mean of the level for which we want a confidence interval.
# xbarj is the mean of the corresponding predictor values.
# It can be a dummy value if type = 0.
# xbar is the overall mean of the predictor values.
# It can be a dummy value if type = 0.
# s2unbl is the mean residual sum of squares from a
# unblocked analysis of variance of the data (where the specimens
# were allocated to the treatments via a predictor sort).
# s2bl is the mean residual sum of squares from a standard
# blocked analysis of variance of the data (where the specimens
# were allocated to the treaments via a predictor sort).
# s2anocov is the mean residual sum of squares from an
# analysis of covariance of the data. If no anocov
# has been performed, this will be a dummy value (and
# type must be 0).
# b is the anocov estimate of the coefficient of
# the predictor. If no anocov has been performed,
# this will be a dummy value (and type must be 0).
# conflev is the confidence level for the confidence interval.
# For a 95% confidence level, conflev should be set to .95.
z = qnorm(conflev + (1 - conflev)/2)
if (type == 0) {
# ANOVA based confidence intervals
sum = 0
prod = 1
for (i in 1:numfac) {
sum = sum + kvec[i]
prod = prod*kvec[i]
}
denom = prod*numrep/(kvec[idfac])
# ad hoc df for the s estimate
df = prod*numrep - sum + (numfac - 1)
t = qt(conflev + (1 - conflev)/2,df)
s = sqrt((s2bl + (s2unbl - s2bl)/kvec[idfac])/denom)
zlow = ybarj - z*s
zup = ybarj + z*s
tlow = ybarj - t*s
tup = ybarj + t*s
list("zlow = ",zlow,"zup = ",zup,"tlow = ",tlow,"tup = ",tup)
} else {
# ANOCOV based confidence intervals
sum = 0
prod = 1
for (i in 1:numfac) {
sum = sum + kvec[i]
prod = prod*kvec[i]
}
denom = prod*numrep/(kvec[idfac])
# ad hoc df for the s estimate
df = prod*numrep - sum + (numfac - 1) - 1
t = qt(conflev + (1 - conflev)/2,df)
s = sqrt((s2anocov + (s2unbl - s2bl)/kvec[idfac])/denom)
est = ybarj - b*(xbarj - xbar)
zlow = est - z*s
zup = est + z*s
tlow = est - t*s
tup = est + t*s
list("zlow = ",zlow,"zup = ",zup,"tlow = ",tlow,"tup = ",tup)
}
}