Nothing
R/EDcomp.R
In drc: Analysis of Dose-Response Curves
Defines functions
createsifct
"EDcomp" <-
function(object, percVec, percMat = NULL, compMatch = NULL, od = FALSE, vcov. = vcov, reverse = FALSE,
interval = c("none", "delta", "fieller", "fls"), level = ifelse(!(interval == "none"), 0.95, NULL),
reference = c("control", "upper"), type = c("relative", "absolute"),
display = TRUE, pool = TRUE, logBase = NULL, multcomp = FALSE, ...)
{
## Matching the argument 'method'
interval <- match.arg(interval)
reference <- match.arg(reference)
type <- match.arg(type)
if ( (is.null(logBase)) && (interval == "fls") )
{
stop("Argument 'logBase' not specified for interval = 'fls'")
}
# sifct <- createsifct(object$"fct"$"edfct", logBase, identical(interval, "fls"), object$"indexMat", length(coef(object)))
## Checking contain of percVec vector ... should be numbers between 0 and 100
if ( (type == "relative") && any(percVec<=0 | percVec>=100) )
{
stop("Percentages outside the interval [0, 100] not allowed")
}
if (missing(compMatch)) {matchNames <- FALSE} else {matchNames <- TRUE}
lenPV <- length(percVec)
## Retrieving relevant quantities
indexMat <- object$"indexMat"
# lenEB <- ncol(indexMat)
parmMat <- object$"parmMat"
# strParm <- colnames(parmMat)
# varMat <- vcov(object, od = od, pool = pool)
## compNames <- as.character(strParm) # converting a factor
curveNames <- colnames(object$"parmMat")
options(warn = -1) # switching off warnings caused by coercion in the if statement
if (any(is.na(as.numeric(curveNames))))
{
curveOrder <- order(curveNames)
} else { # if names are numbers then skip re-ordering
curveOrder <- 1:length(curveNames)
}
options(warn = 0) # normalizing behaviour of warnings
strParm0 <- curveNames[curveOrder]
indexMat <- indexMat[, curveOrder, drop = FALSE]
lenEB <- ncol(indexMat)
sifct <- createsifct(object$"fct"$"edfct", logBase, identical(interval, "fls"), indexMat, length(coef(object)))
parmMat <- parmMat[, curveOrder, drop = FALSE]
strParm <- strParm0
# varMat <- vcov(object, od = od, pool = pool)
varMat <- vcov.(object)
## Calculating SI values
numComp <- (lenPV*(lenPV-1)/2)*(lenEB * (lenEB - 1) / 2)
# if (!identical(interval, "none"))
# {
# siMat <- matrix(0, numComp, 3)
# cNames <- c("Estimate", "Lower", "Upper")
#
# } else {
# siMat <- matrix(0, numComp, 4)
# cNames <- c("Estimate", "Std. Error", "t-value", "p-value")
# }
matchVec <- rep(TRUE, numComp)
rNames <- rep("", numComp)
oriMat <- matrix(0, numComp, 2)
degfree <- df.residual(object)
rowIndex <- 1
# require(gtools, quietly = TRUE)
pairsMat <- combinations(lenEB, 2) # canonical "2" as pairs are considered
if (is.null(percMat))
{
percMat <- combinations(lenPV, 2) # canonical "2" as pairs are considered
}
if (reverse)
{
pairsMat <- pairsMat[, 2:1, drop = FALSE]
percMat <- percMat[, 2:1, drop = FALSE]
# strParm <- rev(strParm)
}
# print(strParm)
# print(pairsMat)
appFct1 <- function(percVal)
{
apply(pairsMat, 1, siInner, pVec = percVec[percVal], compMatch = compMatch, object = object, indexMat = indexMat, parmMat = parmMat,
varMat = varMat, level = level, reference = reference, type = type, sifct = sifct, interval = interval, degfree = degfree, logBase = logBase)
}
SImat0 <- matrix(apply(percMat, 1, appFct1), nrow = nrow(pairsMat) * nrow(percMat), byrow = TRUE)
SImat <- SImat0[, 1:4, drop = FALSE]
# print(SImat)
dSImat <- SImat0[, 5:ncol(SImat0), drop = FALSE]
# matchVec[rowIndex] <- (is.null(compMatch) || all(c(strParm[j], strParm[k]) %in% compMatch))
# strParm0 <- sort(colnames(object$"parmMat"))
appFct2 <- function(percVal)
{
apply(pairsMat, 1,
function(indPair, percVal)
{
# paste(strParm0[indPair[1]], "/", strParm0[indPair[2]], ":", percVec[percVal[1]], "/", percVec[percVal[2]], sep = "")
paste(strParm[indPair[1]], "/", strParm[indPair[2]], ":", percVec[percVal[1]], "/", percVec[percVal[2]], sep = "")
}, percVal = percVal)
}
rownames(SImat) <- apply(percMat, 1, appFct2)
appFct3 <- function(percVal)
{
apply(pairsMat, 1,
function(indPair, percVal)
{
# (is.null(compMatch) || all(c(strParm0[indPair[1]], strParm0[indPair[2]]) %in% compMatch))
(is.null(compMatch) || all(c(strParm[indPair[1]], strParm[indPair[2]]) %in% compMatch))
})
}
SImat <- SImat[as.vector(apply(percMat, 1, appFct3)), , drop = FALSE]
if (!identical(interval, "none"))
{
# siMat <- matrix(0, numComp, 3)
SImat <- SImat[, -4, drop = FALSE]
cNames <- c("Estimate", "Lower", "Upper")
} else {
# siMat <- matrix(0, numComp, 4)
cNames <- c("Estimate", "Std. Error", "t-value", "p-value")
}
colnames(SImat) <- cNames
# print(SImat)
ciLabel <- switch(interval,
"delta" = "Delta method",
"tfls" = "To and from log scale",
"fls" = "From log scale",
"fieller" = "Fieller")
# resPrint(SImat, "Estimated ratios of effect doses\n", interval, ciLabel, display = display)
resPrint(SImat, "Estimated ratios of effect doses", interval, ciLabel, display = display)
## invisible(SImat)
# ## require(multcomp, quietly = TRUE)
## invisible(list(SImat, SImultcomp = list(EDest = EDmat[, 1], EDvcov = dEDmat %*% vcMat %*% t(dEDmat))))
## these lines are older
# invisible(list(SIdisplay = SImat, SImultcomp = list(SImat[, 1], dSImat %*% varMat %*% t(dSImat))))
if(multcomp)
{
SImat1 <- SImat[, 1]
namesVec <- names(SImat1)
SImat1VC <- dSImat %*% varMat %*% t(dSImat)
colnames(SImat1VC) <- namesVec
rownames(SImat1VC) <- namesVec
# invisible(list(multcomp = list(SImat[, 1], dSImat %*% varMat %*% t(dSImat))))
invisible(list(multcomp = parm(SImat1, SImat1VC)))
} else {
invisible(SImat)
}
#
# if (FALSE)
# {
# for (i in 1:lenPV)
# {
# for (ii in 1:lenPV)
# {
# if (i>=ii) {next}
# pVec <- percVec[c(i, ii)]
#
# for (j in 1:lenEB)
# {
# for (k in 1:lenEB)
# {
# if (j>=k) {next}
# matchVec[rowIndex] <- (is.null(compMatch) || all(c(strParm[j], strParm[k]) %in% compMatch))
#
# jInd <- j
# kInd <- k
# if (reverse)
# {
# jInd <- k; kInd <- j; pVec <- pVec[c(2, 1)]
# }
#
# parmInd1 <- indexMat[, jInd]
# parmInd2 <- indexMat[, kInd]
#
# splInd <- splitInd(parmInd1, parmInd2)
#
# parmChosen1 <- parmMat[, jInd]
# parmChosen2 <- parmMat[, kInd]
#
# SIeval <-
# sifct(parmChosen1, parmChosen2, pVec,
# splInd[[1]][, 1], splInd[[2]][, 1], splInd[[3]][, 1], splInd[[3]][, 2], reference, type, jInd, kInd, ...)
#
# indInOrder <- c(splInd[[1]][, 2], splInd[[2]][, 2], splInd[[3]][, 3])
#
# SIval <- SIeval$"val" # SIeval[[1]]
# dSIval <- SIeval$"der" # SIeval[[2]]
## print(dSIval)
#
# oriMat[rowIndex, 1] <- SIval
## oriMat[rowIndex, 2] <- sqrt(dSIval %*% varMat[indInOrder, indInOrder] %*% dSIval) # sqrt(dSIval%*%varCov%*%dSIval)
# oriMat[rowIndex, 2] <- sqrt(t(dSIval) %*% varMat %*% dSIval)
#
# siMat[rowIndex, 1] <- SIval
# rNames[rowIndex] <- paste(strParm[jInd], "/", strParm[kInd], ":", pVec[1], "/", pVec[2], sep="")
#
# ## Using t-distribution for continuous data
# ## only under the normality assumption
# if (identical(object$"type", "continuous"))
# {
# qFct <- function(x) {qt(x, degfree)}
# pFct <- function(x) {pt(x, degfree)}
# } else {
# qFct <- qnorm
# pFct <- pnorm
# }
#
# if (identical(interval, "none"))
# {
# siMat[rowIndex, 2] <- oriMat[rowIndex, 2] # sqrt(dSIval%*%varCov%*%dSIval)
#
# ## Testing SI equal to 1
# tempStat <- (siMat[rowIndex, 1] - 1)/siMat[rowIndex, 2]
# siMat[rowIndex, 3] <- tempStat
# siMat[rowIndex, 4] <- pFct(-abs(tempStat)) + (1 - pFct(abs(tempStat)))
# }
# if ( (identical(interval, "delta")) || (identical(interval, "fls")) )
# {
# stErr <- oriMat[rowIndex, 2] # sqrt(derEval%*%varCov%*%derEval)
# tquan <- qFct(1 - (1 - level)/2)
#
# siMat[rowIndex, 2] <- siMat[rowIndex, 1] - tquan * stErr
# siMat[rowIndex, 3] <- siMat[rowIndex, 1] + tquan * stErr
# ciLabel <- "Delta method"
# }
# if (identical(interval, "tfls"))
# {
# lsVal <- log(oriMat[rowIndex, 1])
# lsdVal <- oriMat[rowIndex, 2]/oriMat[rowIndex, 1]
# tquan <- qFct(1 - (1 - level)/2)
#
# siMat[rowIndex, 2] <- exp(lsVal - tquan * lsdVal)
# siMat[rowIndex, 3] <- exp(lsVal + tquan * lsdVal)
# ciLabel <- "To and from log scale"
# }
# if ((!is.null(logBase)) && (identical(interval, "fls")))
# {
# siMat[rowIndex, 1] <- logBase^(siMat[rowIndex, 1])
# siMat[rowIndex, 2] <- logBase^(siMat[rowIndex, 2])
# siMat[rowIndex, 3] <- logBase^(siMat[rowIndex, 3])
# ciLabel <- "From log scale"
# }
# if (identical(interval, "fieller")) # using t-distribution
# {
# vcMat <- matrix(NA, 2, 2)
# vcMat[1, 1] <- SIeval$"der1"%*%varMat[parmInd1, parmInd1]%*%SIeval$"der1"
# vcMat[2, 2] <- SIeval$"der2"%*%varMat[parmInd2, parmInd2]%*%SIeval$"der2"
# vcMat[1, 2] <- SIeval$"der1"%*%varMat[parmInd1, parmInd2]%*%SIeval$"der2"
# vcMat[2, 1] <- vcMat[1, 2]
# muVec <- c(SIeval$"valnum", SIeval$"valden")
#
# siMat[rowIndex, 2:3] <- fieller(muVec, degfree, vcMat, level = level)
# ciLabel <- "Fieller"
# }
#
#
# rowIndex <- rowIndex+1
# }
# }
# }
# }
# dimnames(siMat) <- list(rNames, cNames)
# siMat <- siMat[matchVec, , drop = FALSE]
#
# resPrint(siMat, "Estimated ratios of effect doses\n", interval, ciLabel, display = display)
# if (display)
# {
# cat("\n")
# cat("Estimated ratios of effect doses\n")
# if (!(ci == "none"))
# {
# ciText <- paste("(", ciLabel, "-based confidence interval(s))\n", sep = "")
# cat(ciText)
# }
# cat("\n")
# printCoefmat(siMat)
# }
# invisible(siMat)
}
#SI <- EDcomp
"fieller" <-
function(mu, df, vcMat, level = 0.95, finney = FALSE, resVar)
{
tper <- qt(1-0.5*(1-level), df)^2
if (!finney)
{
## Based on the entry on Fieller's theorem
## in Encyclopedia of Statistical Sciences Vol. 3 (1983), p. 86
## essentially same formula as in Finney (see below)
mup <- prod(mu)
fVec0 <- mup - tper*vcMat[1,2]
y2 <- mu[2]^2
fVec <- (fVec0)^2 - (mu[1]^2 - tper*vcMat[1,1])*(y2 - tper*vcMat[2,2])
denom <- y2 - tper*vcMat[2,2]
lowerL <- (fVec0 - sqrt(fVec))/denom
upperL <- (fVec0 + sqrt(fVec))/denom
} else {
## Using the formula
## in Finney: Statistical Method in Biological Assay p. 81 (3rd edition, 1978)
## OOPS: uses the estimated residual variance
fac <- sqrt(tper)*sqrt(resVar)/mu[2]
g <- tper*vcMat[2,2]/(mu[2]^2)
if (g >= 1) {stop("Fieller's theorem not useful!")}
ratio <- mu[1]/mu[2]
v11 <- vcMat[1,1]/(resVar)
v12 <- vcMat[1,2]/(resVar)
v22 <- vcMat[2,2]/(resVar)
innerBr <- g*(v11 - (v12^2)/v22)
inBr <- v11 - 2*ratio*v12 + (ratio^2)*v22 - innerBr
firstTerm <- ratio - g*vcMat[1,2]/vcMat[2,2]
secondTerm <- fac*sqrt(inBr)
denom <- 1 - g
lowerL <- (firstTerm - secondTerm)/denom
upperL <- (firstTerm + secondTerm)/denom
}
return(c(lowerL, upperL))
}
"splitInd" <- function(ind1, ind2)
{
matchVec1 <- ind1%in%ind2
matchVec2 <- ind2%in%ind1
# inCommon <- list(pos1 = (1:length(ind1))[matchVec1], pos2 = (1:length(ind2))[matchVec2], val = ind1[matchVec1])
#
lmv1 <- sum(matchVec1)
if (lmv1 > 0.01)
{
inCommon <- matrix( c( (1:length(ind1))[matchVec1], (1:length(ind2))[matchVec2], ind1[matchVec1]), lmv1, 3)
} else {
inCommon <- NULL
}
# only1 <- list(pos = (1:length(ind1))[!matchVec1], val = ind1[!matchVec1])
only1 <- matrix( c( (1:length(ind1))[!matchVec1], ind1[!matchVec1] ), sum(!matchVec1), 2)
# only2 <- list(pos = (1:length(ind2))[!matchVec2], val = ind2[!matchVec2])
only2 <- matrix( c( (1:length(ind2))[!matchVec2], ind2[!matchVec2] ), sum(!matchVec2), 2)
return(list(only1, only2, inCommon))
}
createsifct <- function(edfct, logBase = NULL, fls = FALSE, indexMat, lenCoef)
{
if (is.null(edfct))
{
stop("SI values cannot be calculated")
} else {
if (!fls)
{
if (is.null(logBase)) # this clause has been updated October 12 2010
{
"sifct" <- function(parm1, parm2, pair, jInd, kInd, reference, type, ...)
## ind1, ind2, cmonInd1, cmonInd2 not used
{
# print(parm1)
# print(parm2)
# print(indexMat)
ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
ED1v <- ED1[[1]]
ED1d <- rep(0, lenCoef)
# print(indexMat[, jInd])
# print(ED1[[2]])
ED1d[indexMat[, jInd]] <- ED1[[2]]
# print(ED1v)
# print(ED1d)
ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
ED2v <- ED2[[1]]
ED2d <- rep(0, lenCoef)
ED2d[indexMat[, kInd]] <- ED2[[2]]
# print(ED2v)
# print(ED2d)
SIpair <- ED1v / ED2v # calculating the SI value
SIder <- (ED1d - SIpair * ED2d) / ED2v # calculating the derivative of SI
# SIder1 <- ED1d/ED2v
# SIder2 <- (-ED2d/ED2v)*SIpair
# SIder12 <- commonParm(SIder1, SIder2, cmonInd1, cmonInd2)
# SIder12 <- ED1d/ED2v - (ED2d/ED2v)*SIpair
# return(list(val = SIpair, der = c(SIder1[ind1], SIder2[ind2], SIder12),
return(list(val = SIpair, der = SIder,
der1 = ED1d, der2 = ED2d, valnum = ED1v, valden = ED2v))
}
} else {
"sifct" <- function(parm1, parm2, pair, jInd, kInd, reference, type, ...)
{
ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
ED1v <- ED1[[1]]
ED1d <- rep(0, lenCoef)
ED1d[indexMat[, jInd]] <- ED1[[2]]
ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
ED2v <- ED2[[1]]
ED2d <- rep(0, lenCoef)
ED2d[indexMat[, kInd]] <- ED2[[2]]
# ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
# ED1v <- ED1[[1]]
# ED1d <- ED1[[2]]
# ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
# ED2v <- ED2[[1]]
# ED2d <- ED2[[2]]
SIpair <- logBase^(ED1v - ED2v) # calculating the SI value
SIder <- SIpair * log(logBase) * (ED1d - ED2d)
# SIder1 <- SIpair*log(logBase)*ED1d
# SIder2 <- SIpair*log(logBase)*(-ED2d)
# SIder12 <- commonParm(SIder1, SIder2, cmonInd1, cmonInd2)
# return(list(val = SIpair, der = c(SIder1[ind1], SIder2[ind2], SIder12),
return(list(val = SIpair, der = SIder,
der1 = (log(logBase)*logBase^ED1v)*ED1d, der2 = (log(logBase)*logBase^ED2v)*ED2d,
valnum = logBase^ED1v, valden = logBase^ED2v))
}
}
} else {
"sifct" <- function(parm1, parm2, pair, jInd, kInd, reference, type, ...)
{
ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
ED1v <- ED1[[1]]
ED1d <- rep(0, lenCoef)
ED1d[indexMat[, jInd]] <- ED1[[2]]
ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
ED2v <- ED2[[1]]
ED2d <- rep(0, lenCoef)
ED2d[indexMat[, kInd]] <- ED2[[2]]
# ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
# ED1v <- ED1[[1]]
# ED1d <- ED1[[2]]
# ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
# ED2v <- ED2[[1]]
# ED2d <- ED2[[2]]
SIpair <- ED1v - ED2v # calculating the log SI value
SIder <- ED1d - ED2d
# SIder1 <- ED1d
# SIder2 <- -ED2d
# SIder12 <- commonParm(SIder1, SIder2, cmonInd1, cmonInd2)
# return(list(val = SIpair, der = c(SIder1[ind1], SIder2[ind2], SIder12),
return(list(val = SIpair, der = SIder,
der1 = ED1d, der2 = ED2d, valnum = ED1v, valden = ED2v))
}
}
return(sifct)
}
}
#commonParm <- function(SIder1, SIder2, cmonInd1, cmonInd2)
#{
# lind1 <- length(cmonInd1)
# retVec <- rep(NA, lind1)
# for (i in 1:lind1)
# {
# retVec[i] <- SIder1[cmonInd1[i]] + SIder2[cmonInd2[i]]
#
# }
# return(retVec)
#}
#
Try the drc package in your browser
Any scripts or data that you put into this service are public.
drc documentation built on May 1, 2019, 8:43 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions createsifct
"EDcomp" <-
function(object, percVec, percMat = NULL, compMatch = NULL, od = FALSE, vcov. = vcov, reverse = FALSE,
interval = c("none", "delta", "fieller", "fls"), level = ifelse(!(interval == "none"), 0.95, NULL),
reference = c("control", "upper"), type = c("relative", "absolute"),
display = TRUE, pool = TRUE, logBase = NULL, multcomp = FALSE, ...)
{
## Matching the argument 'method'
interval <- match.arg(interval)
reference <- match.arg(reference)
type <- match.arg(type)
if ( (is.null(logBase)) && (interval == "fls") )
{
stop("Argument 'logBase' not specified for interval = 'fls'")
}
# sifct <- createsifct(object$"fct"$"edfct", logBase, identical(interval, "fls"), object$"indexMat", length(coef(object)))
## Checking contain of percVec vector ... should be numbers between 0 and 100
if ( (type == "relative") && any(percVec<=0 | percVec>=100) )
{
stop("Percentages outside the interval [0, 100] not allowed")
}
if (missing(compMatch)) {matchNames <- FALSE} else {matchNames <- TRUE}
lenPV <- length(percVec)
## Retrieving relevant quantities
indexMat <- object$"indexMat"
# lenEB <- ncol(indexMat)
parmMat <- object$"parmMat"
# strParm <- colnames(parmMat)
# varMat <- vcov(object, od = od, pool = pool)
## compNames <- as.character(strParm) # converting a factor
curveNames <- colnames(object$"parmMat")
options(warn = -1) # switching off warnings caused by coercion in the if statement
if (any(is.na(as.numeric(curveNames))))
{
curveOrder <- order(curveNames)
} else { # if names are numbers then skip re-ordering
curveOrder <- 1:length(curveNames)
}
options(warn = 0) # normalizing behaviour of warnings
strParm0 <- curveNames[curveOrder]
indexMat <- indexMat[, curveOrder, drop = FALSE]
lenEB <- ncol(indexMat)
sifct <- createsifct(object$"fct"$"edfct", logBase, identical(interval, "fls"), indexMat, length(coef(object)))
parmMat <- parmMat[, curveOrder, drop = FALSE]
strParm <- strParm0
# varMat <- vcov(object, od = od, pool = pool)
varMat <- vcov.(object)
## Calculating SI values
numComp <- (lenPV*(lenPV-1)/2)*(lenEB * (lenEB - 1) / 2)
# if (!identical(interval, "none"))
# {
# siMat <- matrix(0, numComp, 3)
# cNames <- c("Estimate", "Lower", "Upper")
#
# } else {
# siMat <- matrix(0, numComp, 4)
# cNames <- c("Estimate", "Std. Error", "t-value", "p-value")
# }
matchVec <- rep(TRUE, numComp)
rNames <- rep("", numComp)
oriMat <- matrix(0, numComp, 2)
degfree <- df.residual(object)
rowIndex <- 1
# require(gtools, quietly = TRUE)
pairsMat <- combinations(lenEB, 2) # canonical "2" as pairs are considered
if (is.null(percMat))
{
percMat <- combinations(lenPV, 2) # canonical "2" as pairs are considered
}
if (reverse)
{
pairsMat <- pairsMat[, 2:1, drop = FALSE]
percMat <- percMat[, 2:1, drop = FALSE]
# strParm <- rev(strParm)
}
# print(strParm)
# print(pairsMat)
appFct1 <- function(percVal)
{
apply(pairsMat, 1, siInner, pVec = percVec[percVal], compMatch = compMatch, object = object, indexMat = indexMat, parmMat = parmMat,
varMat = varMat, level = level, reference = reference, type = type, sifct = sifct, interval = interval, degfree = degfree, logBase = logBase)
}
SImat0 <- matrix(apply(percMat, 1, appFct1), nrow = nrow(pairsMat) * nrow(percMat), byrow = TRUE)
SImat <- SImat0[, 1:4, drop = FALSE]
# print(SImat)
dSImat <- SImat0[, 5:ncol(SImat0), drop = FALSE]
# matchVec[rowIndex] <- (is.null(compMatch) || all(c(strParm[j], strParm[k]) %in% compMatch))
# strParm0 <- sort(colnames(object$"parmMat"))
appFct2 <- function(percVal)
{
apply(pairsMat, 1,
function(indPair, percVal)
{
# paste(strParm0[indPair[1]], "/", strParm0[indPair[2]], ":", percVec[percVal[1]], "/", percVec[percVal[2]], sep = "")
paste(strParm[indPair[1]], "/", strParm[indPair[2]], ":", percVec[percVal[1]], "/", percVec[percVal[2]], sep = "")
}, percVal = percVal)
}
rownames(SImat) <- apply(percMat, 1, appFct2)
appFct3 <- function(percVal)
{
apply(pairsMat, 1,
function(indPair, percVal)
{
# (is.null(compMatch) || all(c(strParm0[indPair[1]], strParm0[indPair[2]]) %in% compMatch))
(is.null(compMatch) || all(c(strParm[indPair[1]], strParm[indPair[2]]) %in% compMatch))
})
}
SImat <- SImat[as.vector(apply(percMat, 1, appFct3)), , drop = FALSE]
if (!identical(interval, "none"))
{
# siMat <- matrix(0, numComp, 3)
SImat <- SImat[, -4, drop = FALSE]
cNames <- c("Estimate", "Lower", "Upper")
} else {
# siMat <- matrix(0, numComp, 4)
cNames <- c("Estimate", "Std. Error", "t-value", "p-value")
}
colnames(SImat) <- cNames
# print(SImat)
ciLabel <- switch(interval,
"delta" = "Delta method",
"tfls" = "To and from log scale",
"fls" = "From log scale",
"fieller" = "Fieller")
# resPrint(SImat, "Estimated ratios of effect doses\n", interval, ciLabel, display = display)
resPrint(SImat, "Estimated ratios of effect doses", interval, ciLabel, display = display)
## invisible(SImat)
# ## require(multcomp, quietly = TRUE)
## invisible(list(SImat, SImultcomp = list(EDest = EDmat[, 1], EDvcov = dEDmat %*% vcMat %*% t(dEDmat))))
## these lines are older
# invisible(list(SIdisplay = SImat, SImultcomp = list(SImat[, 1], dSImat %*% varMat %*% t(dSImat))))
if(multcomp)
{
SImat1 <- SImat[, 1]
namesVec <- names(SImat1)
SImat1VC <- dSImat %*% varMat %*% t(dSImat)
colnames(SImat1VC) <- namesVec
rownames(SImat1VC) <- namesVec
# invisible(list(multcomp = list(SImat[, 1], dSImat %*% varMat %*% t(dSImat))))
invisible(list(multcomp = parm(SImat1, SImat1VC)))
} else {
invisible(SImat)
}
#
# if (FALSE)
# {
# for (i in 1:lenPV)
# {
# for (ii in 1:lenPV)
# {
# if (i>=ii) {next}
# pVec <- percVec[c(i, ii)]
#
# for (j in 1:lenEB)
# {
# for (k in 1:lenEB)
# {
# if (j>=k) {next}
# matchVec[rowIndex] <- (is.null(compMatch) || all(c(strParm[j], strParm[k]) %in% compMatch))
#
# jInd <- j
# kInd <- k
# if (reverse)
# {
# jInd <- k; kInd <- j; pVec <- pVec[c(2, 1)]
# }
#
# parmInd1 <- indexMat[, jInd]
# parmInd2 <- indexMat[, kInd]
#
# splInd <- splitInd(parmInd1, parmInd2)
#
# parmChosen1 <- parmMat[, jInd]
# parmChosen2 <- parmMat[, kInd]
#
# SIeval <-
# sifct(parmChosen1, parmChosen2, pVec,
# splInd[[1]][, 1], splInd[[2]][, 1], splInd[[3]][, 1], splInd[[3]][, 2], reference, type, jInd, kInd, ...)
#
# indInOrder <- c(splInd[[1]][, 2], splInd[[2]][, 2], splInd[[3]][, 3])
#
# SIval <- SIeval$"val" # SIeval[[1]]
# dSIval <- SIeval$"der" # SIeval[[2]]
## print(dSIval)
#
# oriMat[rowIndex, 1] <- SIval
## oriMat[rowIndex, 2] <- sqrt(dSIval %*% varMat[indInOrder, indInOrder] %*% dSIval) # sqrt(dSIval%*%varCov%*%dSIval)
# oriMat[rowIndex, 2] <- sqrt(t(dSIval) %*% varMat %*% dSIval)
#
# siMat[rowIndex, 1] <- SIval
# rNames[rowIndex] <- paste(strParm[jInd], "/", strParm[kInd], ":", pVec[1], "/", pVec[2], sep="")
#
# ## Using t-distribution for continuous data
# ## only under the normality assumption
# if (identical(object$"type", "continuous"))
# {
# qFct <- function(x) {qt(x, degfree)}
# pFct <- function(x) {pt(x, degfree)}
# } else {
# qFct <- qnorm
# pFct <- pnorm
# }
#
# if (identical(interval, "none"))
# {
# siMat[rowIndex, 2] <- oriMat[rowIndex, 2] # sqrt(dSIval%*%varCov%*%dSIval)
#
# ## Testing SI equal to 1
# tempStat <- (siMat[rowIndex, 1] - 1)/siMat[rowIndex, 2]
# siMat[rowIndex, 3] <- tempStat
# siMat[rowIndex, 4] <- pFct(-abs(tempStat)) + (1 - pFct(abs(tempStat)))
# }
# if ( (identical(interval, "delta")) || (identical(interval, "fls")) )
# {
# stErr <- oriMat[rowIndex, 2] # sqrt(derEval%*%varCov%*%derEval)
# tquan <- qFct(1 - (1 - level)/2)
#
# siMat[rowIndex, 2] <- siMat[rowIndex, 1] - tquan * stErr
# siMat[rowIndex, 3] <- siMat[rowIndex, 1] + tquan * stErr
# ciLabel <- "Delta method"
# }
# if (identical(interval, "tfls"))
# {
# lsVal <- log(oriMat[rowIndex, 1])
# lsdVal <- oriMat[rowIndex, 2]/oriMat[rowIndex, 1]
# tquan <- qFct(1 - (1 - level)/2)
#
# siMat[rowIndex, 2] <- exp(lsVal - tquan * lsdVal)
# siMat[rowIndex, 3] <- exp(lsVal + tquan * lsdVal)
# ciLabel <- "To and from log scale"
# }
# if ((!is.null(logBase)) && (identical(interval, "fls")))
# {
# siMat[rowIndex, 1] <- logBase^(siMat[rowIndex, 1])
# siMat[rowIndex, 2] <- logBase^(siMat[rowIndex, 2])
# siMat[rowIndex, 3] <- logBase^(siMat[rowIndex, 3])
# ciLabel <- "From log scale"
# }
# if (identical(interval, "fieller")) # using t-distribution
# {
# vcMat <- matrix(NA, 2, 2)
# vcMat[1, 1] <- SIeval$"der1"%*%varMat[parmInd1, parmInd1]%*%SIeval$"der1"
# vcMat[2, 2] <- SIeval$"der2"%*%varMat[parmInd2, parmInd2]%*%SIeval$"der2"
# vcMat[1, 2] <- SIeval$"der1"%*%varMat[parmInd1, parmInd2]%*%SIeval$"der2"
# vcMat[2, 1] <- vcMat[1, 2]
# muVec <- c(SIeval$"valnum", SIeval$"valden")
#
# siMat[rowIndex, 2:3] <- fieller(muVec, degfree, vcMat, level = level)
# ciLabel <- "Fieller"
# }
#
#
# rowIndex <- rowIndex+1
# }
# }
# }
# }
# dimnames(siMat) <- list(rNames, cNames)
# siMat <- siMat[matchVec, , drop = FALSE]
#
# resPrint(siMat, "Estimated ratios of effect doses\n", interval, ciLabel, display = display)
# if (display)
# {
# cat("\n")
# cat("Estimated ratios of effect doses\n")
# if (!(ci == "none"))
# {
# ciText <- paste("(", ciLabel, "-based confidence interval(s))\n", sep = "")
# cat(ciText)
# }
# cat("\n")
# printCoefmat(siMat)
# }
# invisible(siMat)
}
#SI <- EDcomp
"fieller" <-
function(mu, df, vcMat, level = 0.95, finney = FALSE, resVar)
{
tper <- qt(1-0.5*(1-level), df)^2
if (!finney)
{
## Based on the entry on Fieller's theorem
## in Encyclopedia of Statistical Sciences Vol. 3 (1983), p. 86
## essentially same formula as in Finney (see below)
mup <- prod(mu)
fVec0 <- mup - tper*vcMat[1,2]
y2 <- mu[2]^2
fVec <- (fVec0)^2 - (mu[1]^2 - tper*vcMat[1,1])*(y2 - tper*vcMat[2,2])
denom <- y2 - tper*vcMat[2,2]
lowerL <- (fVec0 - sqrt(fVec))/denom
upperL <- (fVec0 + sqrt(fVec))/denom
} else {
## Using the formula
## in Finney: Statistical Method in Biological Assay p. 81 (3rd edition, 1978)
## OOPS: uses the estimated residual variance
fac <- sqrt(tper)*sqrt(resVar)/mu[2]
g <- tper*vcMat[2,2]/(mu[2]^2)
if (g >= 1) {stop("Fieller's theorem not useful!")}
ratio <- mu[1]/mu[2]
v11 <- vcMat[1,1]/(resVar)
v12 <- vcMat[1,2]/(resVar)
v22 <- vcMat[2,2]/(resVar)
innerBr <- g*(v11 - (v12^2)/v22)
inBr <- v11 - 2*ratio*v12 + (ratio^2)*v22 - innerBr
firstTerm <- ratio - g*vcMat[1,2]/vcMat[2,2]
secondTerm <- fac*sqrt(inBr)
denom <- 1 - g
lowerL <- (firstTerm - secondTerm)/denom
upperL <- (firstTerm + secondTerm)/denom
}
return(c(lowerL, upperL))
}
"splitInd" <- function(ind1, ind2)
{
matchVec1 <- ind1%in%ind2
matchVec2 <- ind2%in%ind1
# inCommon <- list(pos1 = (1:length(ind1))[matchVec1], pos2 = (1:length(ind2))[matchVec2], val = ind1[matchVec1])
#
lmv1 <- sum(matchVec1)
if (lmv1 > 0.01)
{
inCommon <- matrix( c( (1:length(ind1))[matchVec1], (1:length(ind2))[matchVec2], ind1[matchVec1]), lmv1, 3)
} else {
inCommon <- NULL
}
# only1 <- list(pos = (1:length(ind1))[!matchVec1], val = ind1[!matchVec1])
only1 <- matrix( c( (1:length(ind1))[!matchVec1], ind1[!matchVec1] ), sum(!matchVec1), 2)
# only2 <- list(pos = (1:length(ind2))[!matchVec2], val = ind2[!matchVec2])
only2 <- matrix( c( (1:length(ind2))[!matchVec2], ind2[!matchVec2] ), sum(!matchVec2), 2)
return(list(only1, only2, inCommon))
}
createsifct <- function(edfct, logBase = NULL, fls = FALSE, indexMat, lenCoef)
{
if (is.null(edfct))
{
stop("SI values cannot be calculated")
} else {
if (!fls)
{
if (is.null(logBase)) # this clause has been updated October 12 2010
{
"sifct" <- function(parm1, parm2, pair, jInd, kInd, reference, type, ...)
## ind1, ind2, cmonInd1, cmonInd2 not used
{
# print(parm1)
# print(parm2)
# print(indexMat)
ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
ED1v <- ED1[[1]]
ED1d <- rep(0, lenCoef)
# print(indexMat[, jInd])
# print(ED1[[2]])
ED1d[indexMat[, jInd]] <- ED1[[2]]
# print(ED1v)
# print(ED1d)
ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
ED2v <- ED2[[1]]
ED2d <- rep(0, lenCoef)
ED2d[indexMat[, kInd]] <- ED2[[2]]
# print(ED2v)
# print(ED2d)
SIpair <- ED1v / ED2v # calculating the SI value
SIder <- (ED1d - SIpair * ED2d) / ED2v # calculating the derivative of SI
# SIder1 <- ED1d/ED2v
# SIder2 <- (-ED2d/ED2v)*SIpair
# SIder12 <- commonParm(SIder1, SIder2, cmonInd1, cmonInd2)
# SIder12 <- ED1d/ED2v - (ED2d/ED2v)*SIpair
# return(list(val = SIpair, der = c(SIder1[ind1], SIder2[ind2], SIder12),
return(list(val = SIpair, der = SIder,
der1 = ED1d, der2 = ED2d, valnum = ED1v, valden = ED2v))
}
} else {
"sifct" <- function(parm1, parm2, pair, jInd, kInd, reference, type, ...)
{
ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
ED1v <- ED1[[1]]
ED1d <- rep(0, lenCoef)
ED1d[indexMat[, jInd]] <- ED1[[2]]
ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
ED2v <- ED2[[1]]
ED2d <- rep(0, lenCoef)
ED2d[indexMat[, kInd]] <- ED2[[2]]
# ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
# ED1v <- ED1[[1]]
# ED1d <- ED1[[2]]
# ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
# ED2v <- ED2[[1]]
# ED2d <- ED2[[2]]
SIpair <- logBase^(ED1v - ED2v) # calculating the SI value
SIder <- SIpair * log(logBase) * (ED1d - ED2d)
# SIder1 <- SIpair*log(logBase)*ED1d
# SIder2 <- SIpair*log(logBase)*(-ED2d)
# SIder12 <- commonParm(SIder1, SIder2, cmonInd1, cmonInd2)
# return(list(val = SIpair, der = c(SIder1[ind1], SIder2[ind2], SIder12),
return(list(val = SIpair, der = SIder,
der1 = (log(logBase)*logBase^ED1v)*ED1d, der2 = (log(logBase)*logBase^ED2v)*ED2d,
valnum = logBase^ED1v, valden = logBase^ED2v))
}
}
} else {
"sifct" <- function(parm1, parm2, pair, jInd, kInd, reference, type, ...)
{
ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
ED1v <- ED1[[1]]
ED1d <- rep(0, lenCoef)
ED1d[indexMat[, jInd]] <- ED1[[2]]
ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
ED2v <- ED2[[1]]
ED2d <- rep(0, lenCoef)
ED2d[indexMat[, kInd]] <- ED2[[2]]
# ED1 <- edfct(parm1, pair[1], reference = reference, type = type, ...)
# ED1v <- ED1[[1]]
# ED1d <- ED1[[2]]
# ED2 <- edfct(parm2, pair[2], reference = reference, type = type, ...)
# ED2v <- ED2[[1]]
# ED2d <- ED2[[2]]
SIpair <- ED1v - ED2v # calculating the log SI value
SIder <- ED1d - ED2d
# SIder1 <- ED1d
# SIder2 <- -ED2d
# SIder12 <- commonParm(SIder1, SIder2, cmonInd1, cmonInd2)
# return(list(val = SIpair, der = c(SIder1[ind1], SIder2[ind2], SIder12),
return(list(val = SIpair, der = SIder,
der1 = ED1d, der2 = ED2d, valnum = ED1v, valden = ED2v))
}
}
return(sifct)
}
}
#commonParm <- function(SIder1, SIder2, cmonInd1, cmonInd2)
#{
# lind1 <- length(cmonInd1)
# retVec <- rep(NA, lind1)
# for (i in 1:lind1)
# {
# retVec[i] <- SIder1[cmonInd1[i]] + SIder2[cmonInd2[i]]
#
# }
# return(retVec)
#}
#
Try the drc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.