Nothing
#'## Compute CIs for prediction ratios
#'### Functions for computing Ratio CIs
FiellerRatioCI <- function(a, b, V, t)
{
theta <- a/b
g <- (t^2)*V[2,2]/b^2
se <- sqrt(V[1,1] - 2*theta*V[1,2] + theta^2 * V[2,2] - g*(V[1,1]-V[1,2]^2/V[2,2]))
LCL <- (1/(1-g))*(theta- g*V[1,2]/V[2,2] - t/b * se)
UCL <- (1/(1-g))*(theta- g*V[1,2]/V[2,2] + t/b * se)
if (!is.na(UCL) & !is.na(LCL) & UCL < LCL) {
UCL <- NA; LCL <- NA
}
return(list(pred.ratio=theta, upper.Confidence.limit = UCL, lower.Confidence.limit = LCL))
}
makeContrastMatrix <- function(preds, indx, pairs.factor, first.level, second.level)
{
X <- model.matrix(~ Combinations-1, data = preds)
## Form the contrast matrix
#Determine the columns of X for each contrast
colsets <- split(colnames(X), preds[indx], lex.order = TRUE)
prediv <- split(preds, preds[indx], lex.order = TRUE)
locns <- lapply(colsets,
function(x, cols = colnames(X))
{ pair <- na.omit(match(x, cols)); names(pair) <- x; return(pair) })
contr <- mapply(function(pred, cols, pairs.factor, first.level, second.level)
{
names(cols) <- gsub("Combinations", "", names(cols), fixed = TRUE)
# k1 <- (1:length(cols))[grepl(first.level, names(cols), fixed = TRUE)]
# k2 <- (1:length(cols))[grepl(second.level, names(cols), fixed = TRUE)]
k1 <- (1:length(cols))[pred[pairs.factor] == first.level]
k2 <- (1:length(cols))[pred[pairs.factor] == second.level]
if (any(c(length(k1),length(k2)) == 0)) #no contrast
cols <- NA
else
{
if (length(k1) != 1)
stop("There is more than one observation with the level ", first.level, " in combination ", names(cols))
else
{
if (length(k2) != 1)
stop("There is more than one observation with the level ", second.level, " in combination ", names(cols))
else
{
cols[1:length(cols)] <- 0
cols[k1] <- 1
cols[k2] <- -1
}
}
}
return(cols)
}, prediv, locns, MoreArgs = list(pairs.factor = pairs.factor,
first.level = first.level, second.level = second.level),
SIMPLIFY = FALSE)
#Substitute the contrast in one column and NA in the other for each contrast
for (k in 1:length(locns))
{
for (j in 1:length(locns[[k]]))
{
X[locns[[k]][j], locns[[k]][1]] <- contr[[k]][j]
if (j != 1)
X[locns[[k]][j], locns[[k]][j]] <- 0
}
}
#Form contrast matrix
X <- X[, unlist(locns)[order(unlist(locns))]]
which.cols <- unlist(lapply(1:ncol(X),
function(k)
{
if (any(is.na(as.vector(X[,k])))) retain <- FALSE
else retain <- !all(X[,k] == 0)
}))
L <- t(X[, which.cols])
rownames(L) <- gsub(paste0("Combinations"), "", rownames(L), fixed = TRUE)
lev <- levels(preds[[pairs.factor]])[1]
rownames(L) <- gsub(paste0(",",lev), "", rownames(L), fixed = TRUE)
rownames(L) <- gsub(paste0(lev,","), "", rownames(L), fixed = TRUE)
return(L)
}
ratioTransform.alldiffs <- function(alldiffs.obj, ratio.factor,
numerator.levels, denominator.levels,
method = "Fieller", alpha = 0.05,
response = NULL, response.title = NULL,
tables = "predictions",
...)
{
#Check that a valid object of class alldiffs
validalldifs <- validAlldiffs(alldiffs.obj)
if (is.character(validalldifs))
stop(validalldifs)
alldiffs.obj <- renameDiffsAttr(alldiffs.obj)
## Check have the vcov component
if (all(is.null(alldiffs.obj$vcov)))
stop("The alldiffs.obj must have a non-null vcov component - use Vmatrix = TRUE when forming alddifs.obj")
meth.options <- c("Fieller")
method.opt <- meth.options[check.arg.values(method, meth.options)]
options <- c("none", "predictions", "backtransforms", "vcov",
"differences", "p.differences", "sed", "LSD", "all")
opt <- options[unlist(lapply(tables, check.arg.values, options=options))]
#Get a t-value (use z if no tdf attribute)
tdf <- attr(alldiffs.obj, which = "tdf")
if (!is.null(tdf))
{
t <- qt(1 - alpha/2, df = tdf)
} else
t <- qnorm(1 - alpha/2)
if (is.null(response))
response <- attr(alldiffs.obj, which = "response")
if (is.null(response.title))
response.title <- attr(alldiffs.obj, which = "response.title")
classify <- attr(alldiffs.obj, which = "classify")
facs <- fac.getinTerm(classify, rmfunction = TRUE)
alldiffs.obj <- renewClassify(alldiffs.obj, newclassify = classify) #make sure in standard order
tmp <- alldiffs.obj$predictions
tmp$Combinations <- dae::fac.combine(as.list(tmp[facs]), combine.levels = TRUE, sep = ",")
if (length(ratio.factor) != 1 | !(ratio.factor %in% facs))
stop("ratio factor must specify a single factor that is in the classify attribute of the alldiffs object")
if (!all(c(numerator.levels, denominator.levels) %in% levels(tmp[[ratio.factor]])))
stop("Not all numerator.levels and denominator.levels are levels in ", ratio.factor, " in the alldiffs object")
sortFactor <- attr(alldiffs.obj, which = "sortFactor")
indx <- setdiff(facs, ratio.factor)
# if (length(indx) == 1)
# {
# ratios.dat <- as.data.frame(ratios.dat, row.names = NULL)
# names(ratios.dat)[1] <- indx
# }
tmp <- split(tmp, tmp[indx], sep = ",", lex.order = TRUE)
Ratios <- lapply(denominator.levels, function(denom.lev,tmp, t)
{
lapply(numerator.levels, function(num.lev, denom.lev, tmp, t)
{
if (!(num.lev %in% alldiffs.obj$predictions[[ratio.factor]]) ||
!(denom.lev %in% alldiffs.obj$predictions[[ratio.factor]]) ||
num.lev == denom.lev)
ratioCIs <-NULL
else
{
ratioCIs <- do.call(rbind, lapply(tmp, function(preds, num.lev, denom.lev, V, t)
{
indx.vals <- (lapply(preds[indx], function(x) as.character(x)[1])) #current levels of indx factors
ratio.pair <- c(indx.vals, rep(NA,3))
names(ratio.pair) <- c(indx, "pred.ratio", "upper.Confidence.limit", "lower.Confidence.limit")
if (all(c(num.lev, denom.lev) %in% preds[[ratio.factor]]))
{
a <- preds[preds[ratio.factor] == num.lev, "predicted.value"]
b <- preds[preds[ratio.factor] == denom.lev, "predicted.value"]
labs <- preds$Combinations
posns <- unlist(lapply(c(num.lev,denom.lev), grep, x =labs, fixed = TRUE))
lab <- labs[posns]
Vpair <- V[lab,lab]
ratio.pair <- as.data.frame(c(indx.vals, (FiellerRatioCI(a, b, Vpair, t = t))))
}
return(ratio.pair)
}, V = alldiffs.obj$vcov, num.lev = num.lev, denom.lev = denom.lev, t = t))
rownames(ratioCIs) <- NULL
names(ratioCIs)[match("pred.ratio", names(ratioCIs))] <- "predicted.value"
ratioCIs$standard.error <- NA
ratioCIs$est.status <- "Estimable"
ratio.names <- names(ratioCIs)
n <- match("standard.error", ratio.names)
k <- match("upper.Confidence.limit", ratio.names)
if (n-k+1 != 3)
stop("Standard.error and confidence limits not together in ratios data frame")
ratio.names[k:n] <- c("standard.error", "upper.Confidence.limit", "lower.Confidence.limit")
ratioCIs <- ratioCIs[ratio.names]
ratioCIs[indx] <- lapply(indx,
function(fac) new.fac <- factor(ratioCIs[[fac]],
levels = levels(alldiffs.obj$predictions[[fac]])))
ratioCIs <- as.predictions.frame(ratioCIs)
class(ratioCIs) <- c("predictions.frame", "data.frame")
attr(ratioCIs, which = "response") <- response
attr(ratioCIs, which = "response.title") <- response.title
#sort the predictions.frame, if it was previously sorted
if (!is.null(sortFactor))
{
if (sortFactor != ratio.factor)
ratioCIs <- sort(ratioCIs, sortFactor = sortFactor, classify = paste(indx, collapse = ":"),
sortOrder = attr(alldiffs.obj, which = "sortOrder"))
}
}
return(ratioCIs)
}, tmp = tmp, denom.lev = denom.lev, t = t)
}, tmp = tmp, t = t)
Ratios <- unlist(Ratios, recursive = FALSE)
nams <- outer(numerator.levels, denominator.levels, paste, sep = ",")
names(Ratios) <- nams
if ("predictions" == opt)
{
cat("\n\n#### Table(s) of ratio predictions and confidence limits for ", response, "\n\n")
print(Ratios)
}
invisible(Ratios)
}
pairdiffsTransform.alldiffs <- function(alldiffs.obj, pairs.factor, first.levels, second.levels,
Vmatrix = FALSE,
error.intervals = "Confidence",
avsed.tolerance = 0.25, accuracy.threshold = NA,
LSDtype = "overall", LSDsupplied = NULL, LSDby = NULL,
LSDstatistic = "mean", LSDaccuracy = "maxAbsDeviation",
response = NULL, response.title = NULL, tables = "all",
pairwise = TRUE, alpha = 0.05, ...)
{
#Check for deprecated argument meanLSD.type and warn
tempcall <- list(...)
if (length(tempcall))
if ("meanLSD.type" %in% names(tempcall))
stop("meanLSD.type has been deprecated - use LSDtype")
#Check that a valid object of class alldiffs
validalldifs <- validAlldiffs(alldiffs.obj)
if (is.character(validalldifs))
stop(validalldifs)
alldiffs.obj <- renameDiffsAttr(alldiffs.obj)
attr(alldiffs.obj, which = "alpha") <- alpha
#get transform attributes from backtransforms
transform.power <- 1; offset <- 0; scale <- 1; transform.function <- "identity"
if (!is.null(alldiffs.obj$backtransforms))
{
transform.power <- attr(alldiffs.obj$backtransforms, which = "transform.power")
offset <- attr(alldiffs.obj$backtransforms, which = "offset")
scale <- attr(alldiffs.obj$backtransforms, which = "scale")
transform.function <- attr(alldiffs.obj$backtransforms, which = "transform.function")
if (is.null(transform.function))
transform.function <- identity
}
sortFactor <- attr(alldiffs.obj, which = "sortFactor")
int.options <- c("none", "Confidence", "StandardError", "halfLeastSignificant")
int.opt <- int.options[check.arg.values(error.intervals, int.options)]
options <- c("none", "predictions", "backtransforms", "vcov",
"differences", "p.differences", "sed", "LSD", "all")
opt <- options[unlist(lapply(tables, check.arg.values, options=options))]
#Get a t-value (use z if no tdf attribute)
tdf <- attr(alldiffs.obj, which = "tdf")
if (!is.null(tdf))
{
t <- qt(1 - alpha/2, df = tdf)
} else
t <- qnorm(1 - alpha/2)
if (is.null(response))
response <- attr(alldiffs.obj, which = "response")
if (is.null(response.title))
response.title <- attr(alldiffs.obj, which = "response.title")
classify <- attr(alldiffs.obj, which = "classify")
facs <- fac.getinTerm(classify, rmfunction = TRUE)
alldiffs.obj <- renewClassify(alldiffs.obj, newclassify = classify) #make sure in standard order
tmp <- alldiffs.obj$predictions
tmp$Combinations <- dae::fac.combine(as.list(tmp[facs]), combine.levels = TRUE, sep = ",")
if (length(pairs.factor) != 1 | !(pairs.factor %in% facs))
stop("pairs factor must specify a single factor that is in the classify attribute of the alldiffs object")
if (!all(c(first.levels, second.levels) %in% levels(tmp[[pairs.factor]])))
stop("Not all first.levels and second.levels are levels in ", pairs.factor, " in the all diffs object")
indx <- setdiff(facs, pairs.factor)
Diffs <- lapply(second.levels, function(second.lev,tmp)
{
lapply(first.levels, function(first.lev, second.lev, tmp)
{
if (!(first.lev %in% alldiffs.obj$predictions[[pairs.factor]]) ||
!(second.lev %in% alldiffs.obj$predictions[[pairs.factor]]))
diffs <-NULL
else
{
L <- makeContrastMatrix(tmp, indx = indx, pairs.factor = pairs.factor, first.lev, second.lev)
#form the factors indexing the pair-difference predictions
new.facs <- lapply(indx, function(fac, data)
{
levels(tmp[[fac]])
}, data = tmp)
names(new.facs) <- indx
if (length(indx) == 1)
{
pairs.dat <- as.data.frame(rownames(L), stringsAsFactors = FALSE)
names(pairs.dat) <- indx
pairs.dat[indx] <- factor(pairs.dat[[indx]], levels = rownames(L))
}
else
pairs.dat <- dae::fac.uncombine(rownames(L), new.factors = new.facs, sep = ",")
#Calculate the differences for the current pair
diffs <- linTransform(alldiffs.obj, classify = classify,
linear.transformation = L, Vmatrix = Vmatrix,
error.intervals = "Confidence",
avsed.tolerance = avsed.tolerance, accuracy.threshold = accuracy.threshold,
response = response, response.title = response.title,
pairwise = pairwise, alpha = alpha,
tables = "none", ...)
diffs$predictions <- cbind(pairs.dat, diffs$predictions)
diffs$predictions <- diffs$predictions[, -match("Combination", names(diffs$predictions))]
diffs <- renewClassify(diffs, newclassify = paste(indx, collapse = ":"))
diffs <- redoErrorIntervals(diffs, error.intervals = int.opt, alpha = alpha,
avsed.tolerance = avsed.tolerance, accuracy.threshold = accuracy.threshold,
LSDtype = LSDtype, LSDby = LSDby, LSDsupplied = LSDsupplied,
LSDstatistic = LSDstatistic, LSDaccuracy = LSDaccuracy)
# if (!is.null(diffs$backtransforms))
# {
# diffs$backtransforms <- cbind(pairs.dat, diffs$backtransforms)
# diffs$backtransforms <- diffs$backtransforms[, -match("Combination", names(diffs$backtransforms))]
# #Find missing attributes in new alldiffs.obj and add them back in
# newattr <- attributes(diffs$backtransforms)
# back.attr <- attributes(alldiffs.obj$backtransforms)
# back.attr <- back.attr[names(back.attr)[!(names(back.attr) %in% names(newattr))]]
# if (length(back.attr) > 0)
# {
# newattr <- c(newattr,back.attr)
# attributes(diffs$backtransforms) <- newattr
# }
# }
#sort the alldiffs if it was previously sorted
if (!is.null(sortFactor))
{
if (sortFactor != pairs.factor)
diffs <- sort(diffs, sortFactor = sortFactor,
sortOrder = attr(diffs, which = "sortOrder"))
}
if (opt != "none")
print(diffs, which = opt)
return(diffs)
}
}, tmp = tmp, second.lev = second.lev)
}, tmp = tmp)
Diffs <- unlist(Diffs, recursive = FALSE)
nams <- outer(first.levels, second.levels, paste, sep = ",")
names(Diffs) <- nams
invisible(Diffs)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.