summary.braidrm <- function(object, ...) {
if (is.null(object$ciPass) || !object$ciPass) {
TAB <- coef(object)
} else {
cv <- object$ciVec
TAB <- cbind(CILow = cv[seq(1,length(cv),by=2)],Estimate = coef(object),CIHigh = cv[seq(2,length(cv),by=2)])
}
res <- list(call=object$call,coefficients=TAB)
class(res) <- "summary.braidrm"
return(res)
}
print.summary.braidrm <- function(x, ...) {
cat("Call:\n")
print(x$call)
cat("\n")
print(round(x$coefficients,digits=4))
}
print.braidrm <- function(x, ...) {
cat("Call:\n")
print(x$call)
cat("\nCoefficients:\n")
print(x$fullpar)
}
braidrm <- function(model,data,getCIs=TRUE,fixed="kappa2",startparv=NULL,llims=NULL,ulims=NULL,...) UseMethod("braidrm")
braidrm.default <- function(model,data,getCIs=TRUE,fixed="kappa2",startparv=NULL,llims=NULL,ulims=NULL,...) {
concs <- model
act <- data
if (ncol(concs)!=2) { stop("Parameter 'concs' must be an array with two columns.") }
conc1 <- as.vector(concs[,1])
conc2 <- as.vector(concs[,2])
act <- as.vector(act)
if (length(conc1)!=length(conc2) | length(conc1)!=length(act)) {
stop("Parameters 'concs', and 'act' must contain the same number of elements.")
}
cdef <- startparv
cfixed <- fixed
for (iter in 1:5) {
nls <- try(fitBRAIDrsm(conc1,conc2,act,fixed=cfixed,def=cdef,llims=llims,ulims=ulims),silent=TRUE)
if (class(nls) != "try-error" && nls$convergence==0) { break }
else if (iter==5) { stop(paste("Convergence error:",nls$message)) }
cdef <- nls$fullpar
cdef[1:5] <- log(cdef[1:5])
cdef[1:7] <- cdef[1:7]+runif(7,min=-0.01,max=0.01)
# cdef[8:10] <- mean(cdef[8:10])+runif(1,min=-0.01,max=0.01)
cdef[1:5] <- exp(cdef[1:5])
cfixed <- nls$fixed
if ((is.character(fixed)||length(fixed)<10) && !is.na(cfixed[10])) { cfixed[10] <- cdef[10] }
cdef <- cdef[which(is.na(cfixed))]
cdef <- pmax(cdef,nls$mlims[1,])
cdef <- pmin(cdef,nls$mlims[2,])
}
ostart <- nls$odef
nfixed <- nls$fixed
fullpar <- nls$fullpar
nfixed[which(!is.na(nfixed))] <- fullpar[which(!is.na(nfixed))]
ccoef <- fullpar[which(is.na(nfixed))]
pred_act <- evalBRAIDrsm(conc1,conc2,fullpar)
bfit <- list(conc1=conc1,conc2=conc2,act=act,fitted.values=pred_act,residuals=(pred_act-act),
ostart=ostart,fixed=nfixed,mlims=nls$mlims,coefficients=ccoef,fullpar=fullpar,
convergence=nls$convergence,message=nls$message,call=match.call())
cfnames <- c("IDMA","IDMB","na","nb","delta","kappa","E0","EfA","EfB","EfAB")
names(bfit$coefficients) <- cfnames[which(is.na(nfixed))]
class(bfit) <- "braidrm"
if (getCIs) { bfit <- getBRAIDbootstrap(bfit) }
return(bfit)
}
braidrm.formula <- function(model,data,...) {
mf <- model.frame(formula=model, data=data)
concs <- model.matrix(attr(mf, "terms"), data=mf)
tms <- attr(concs,"assign")
for (i in seq(length(tms),1,by=-1)) {
if (tms[i]==0) { concs <- concs[,-i] }
}
act <- model.response(mf)
bfit <- braidrm.default(concs,act,...)
bfit$call <- match.call()
return(bfit)
}
getBRAIDbootstrap <- function(bfit,ciLevs=c(0.025,0.975),numBoot=NULL) {
if (class(bfit)!="braidrm") { stop("Input 'bfit' must be of class 'braidrm'.") }
if (!is.null(bfit$ciPass)) {
warning("Input already has bootstrapped coefficients (or has failed bootstrapping).")
return(bfit)
}
prcBootPass <- 0.5
if (is.null(numBoot)) { numBoot <- max(min(10/(1-ciLevs[2]+ciLevs[1]),1000),100) }
bCoefs <- array(coef(bfit),dim=c(1,length(coef(bfit))))
for (i in 1:numBoot) {
bact <- fitted(bfit) + sample(resid(bfit),length(resid(bfit)),replace=TRUE)
nls <- try(fitBRAIDrsm(bfit$conc1,bfit$conc2,bact,fixed=bfit$fixed,
def=coef(bfit),llims=bfit$mlims[1,],ulims=bfit$mlims[2,]),silent=TRUE)
if (class(nls)!="try-error" && nls$convergence==0 && !(TRUE %in% (is.infinite(nls$par) | is.nan(nls$par)))) {
cpar <- nls$fullpar[which(is.na(bfit$fixed))]
bCoefs <- rbind(bCoefs,array(cpar,dim=c(1,length(cpar))))
}
}
nbfit <- bfit
if (nrow(bCoefs)<numBoot*prcBootPass) {
nbfit$ciPass <- FALSE
} else {
qMat <- apply(bCoefs,2,quantile,probs=ciLevs)
nbfit <- c(nbfit,list(ciPass=TRUE,ciLevs=ciLevs,bCoefs=bCoefs,ciVec=as.vector(qMat)))
}
class(nbfit) <- "braidrm"
return(nbfit)
}
calcBRAIDconfint <- function(bfit,parfunc,civals=NULL) {
if (class(bfit)!="braidrm") { stop("Input 'bfit' must be of class 'braidrm'.") }
if (class(parfunc)!="function") { stop("Input 'parfunc' must be a function of one variable.") }
if (is.null(bfit$ciPass) || !bfit$ciPass) { stop("Input 'bfit' must have bootstrapped coefficients.") }
if (is.null(civals)) { civals <- bfit$ciLevs }
ostart <- bfit$ostart
fixed <- bfit$fixed
fullpar <- bfit$fullpar
isIncr <- sign(fullpar[10]-fullpar[7])
parv2fullpar <- function(parv) {
tfullpar <- fullpar
tfullpar[which(is.na(fixed))] <- parv
if (is.na(fixed[10])) {
if (is.na(fixed[8])||is.na(fixed[9])) {
if (!is.na(fixed[8]) && ostart[8]==ostart[10]) { tfullpar[8] <- tfullpar[10] }
else if (!is.na(fixed[9]) && ostart[9]==ostart[10]) { tfullpar[9] <- tfullpar[10] }
else { tfullpar[10] <- tfullpar[10]+isIncr*max(isIncr*tfullpar[8:9]) }
}
else if (ostart[8]==ostart[10] && ostart[9]==ostart[10]) { tfullpar[8:9] <- tfullpar[10] }
} else if ((is.na(fixed[8]) || is.na(fixed[9])) && ostart[10]==isIncr*max(isIncr*ostart[8:9])) {
tfullpar[10] <- isIncr*max(isIncr*tfullpar[8:9])
}
return(tfullpar)
}
outval <- parfunc(fullpar)
outmat <- outval
if (length(outmat)==1) { outmat <- array(outmat,dim=c(1,1)) }
for (b in 1:nrow(bfit$bCoefs)) { outmat <- cbind(outmat,parfunc(parv2fullpar(bfit$bCoefs[b,]))) }
outci <- apply(outmat,1,quantile,probs=civals)
fullout <- cbind(as.vector(outci[1,]),outval,as.vector(outci[2,]))
return(fullout)
}
# Params 1 and 2: IDMA and IDMB
# Params 3 and 4: na and nb
# Param 5: delta
# Param 6: kappa
# Param 7: E0
# Params 8 and 9: EA and EB
# Param 10: EAB
evalBRAIDrsm <- function(DA,DB,parv) { return(evalBRAIDrsm_int(DA,DB,parv)) }
evalBRAIDrsm_int <- function(DA,DB,parv,fixed=c(NA,NA,NA,NA,1,NA,NA,NA,NA,NA),calcderivs=FALSE) {
concA <- DA
concB <- DB
if (abs(parv[8]-parv[7])>abs(parv[10]-parv[7]) || abs(parv[9]-parv[7])>abs(parv[10]-parv[7])) {
stop("Invalid parameter values.")
}
if (parv[8]==parv[10]) {
Aden <- 1
Amod <- (concA/parv[1])^(sqrt(parv[3]/parv[4])/parv[5])
} else {
Aden <- 1+(1-((parv[8]-parv[7])/(parv[10]-parv[7])))*((concA/parv[1])^parv[3])
Amod <- ((parv[8]-parv[7])/(parv[10]-parv[7]))*((concA/parv[1])^parv[3])
iinds <- which(is.infinite(Aden)|is.infinite(Amod))
Aden[iinds] <- 1-(parv[8]-parv[7])/(parv[10]-parv[7])
Amod[iinds] <- (parv[8]-parv[7])/(parv[10]-parv[7])
Amod <- (Amod/Aden)^(1/(parv[5]*sqrt(parv[3]*parv[4])))
Aden[iinds] <- Inf
}
if (parv[9]==parv[10]) {
Bden <- 1
Bmod <- (concB/parv[2])^(sqrt(parv[4]/parv[3])/parv[5])
} else {
Bden <- 1+(1-((parv[9]-parv[7])/(parv[10]-parv[7])))*((concB/parv[2])^parv[4])
Bmod <- ((parv[9]-parv[7])/(parv[10]-parv[7]))*((concB/parv[2])^parv[4])
iinds <- which(is.infinite(Bden)|is.infinite(Bmod))
Bden[iinds] <- 1-(parv[9]-parv[7])/(parv[10]-parv[7])
Bmod[iinds] <- (parv[9]-parv[7])/(parv[10]-parv[7])
Bmod <- (Bmod/Bden)^(1/(parv[5]*sqrt(parv[3]*parv[4])))
Bden[iinds] <- Inf
}
ABrad <- sqrt(Amod*Bmod)
if (parv[6]!=0) { ABrad[which(Amod==0|Bmod==0)] <- 0
} else { ABrad[which(Amod==0|Bmod==0|ABrad==Inf)] <- 0 }
Exp <- Amod+Bmod+parv[6]*ABrad
Exp[is.infinite(Amod)|is.infinite(Bmod)] <- Inf
Den <- 1+Exp^(-parv[5]*sqrt(parv[3]*parv[4]))
Ei <- parv[7]+(parv[10]-parv[7])/Den
if (!calcderivs) { return(Ei) }
Emod <- (parv[10]-parv[7])*(Exp^(-parv[5]*sqrt(parv[3]*parv[4])-1))/(Den^2)
Emod[which(Exp==0)] <- 0
Emodz <- which(Emod==0 | is.infinite(Emod))
derivs <- array(0,dim=c(length(Ei),10))
if (is.na(fixed[1])) { derivs[,1] <- -parv[3]*Emod*(Amod+parv[6]*ABrad/2)/(Aden*parv[1]) }
if (is.na(fixed[2])) { derivs[,2] <- -parv[4]*Emod*(Bmod+parv[6]*ABrad/2)/(Bden*parv[2]) }
if (is.na(fixed[3]) || is.na(fixed[4]) || is.na(fixed[5])) {
AlogA <- Amod*log(Amod)
AlogA[which(Amod==0)] <- 0
BlogB <- Bmod*log(Bmod)
BlogB[which(Bmod==0)] <- 0
RlogR <- ABrad*log(ABrad)
RlogR[which(ABrad==0)] <- 0
ElogE <- Exp*log(Exp)
ElogE[which(Exp==0)] <- 0
derivs[,5] <- -sqrt(parv[3]*parv[4])*Emod*(AlogA+BlogB+parv[6]*RlogR-ElogE)
if (is.na(fixed[3])) {
derivs[,3] <- Emod*(Amod+parv[6]*ABrad/2)*log(concA/parv[1])/Aden
derivs[which(concA==0),3] <- 0
derivs[,3] <- derivs[,3]+parv[5]*derivs[,5]/(2*parv[3])
}
if (is.na(fixed[4])) {
derivs[,4] <- Emod*(Bmod+parv[6]*ABrad/2)*log(concB/parv[2])/Bden
derivs[which(concB==0),4] <- 0
derivs[,4] <- derivs[,4]+parv[5]*derivs[,5]/(2*parv[4])
}
}
if (is.na(fixed[6])) { derivs[,6] <- parv[5]*sqrt(parv[3]*parv[4])*Emod*ABrad }
derivs[Emodz,1:6] <- 0
if (is.na(fixed[7]) || is.na(fixed[8]) || is.na(fixed[9]) || is.na(fixed[10])) {
if (is.na(fixed[10]) && !is.na(fixed[8]) && !is.na(fixed[9]) && parv[8]==parv[10] && parv[9]==parv[10]) {
derivs[,10] <- 1/Den
} else if (xor(is.na(fixed[8]),is.na(fixed[10])) && parv[8]==parv[10]) {
derivs[,9] <- Emod*(Bmod+parv[6]*ABrad/2)*(1+(concB/parv[2])^parv[4])/(Bden*(parv[9]-parv[7]))
derivs[Emodz,9] <- 0
if (is.na(fixed[8])) { derivs[,8] <- 1/Den-derivs[,9]*(parv[9]-parv[7])/(parv[8]-parv[7]) }
else { derivs[,10] <- 1/Den-derivs[,9]*(parv[9]-parv[7])/(parv[10]-parv[7]) }
} else if (xor(is.na(fixed[9]),is.na(fixed[10])) && parv[9]==parv[10]) {
derivs[,8] <- Emod*(Amod+parv[6]*ABrad/2)*(1+(concA/parv[1])^parv[3])/(Aden*(parv[8]-parv[7]))
derivs[Emodz,8] <- 0
if (is.na(fixed[9])) { derivs[,9] <- 1/Den-derivs[,8]*(parv[8]-parv[7])/(parv[9]-parv[7]) }
else { derivs[,10] <- 1/Den-derivs[,8]*(parv[8]-parv[7])/(parv[10]-parv[7]) }
} else {
derivs[,8] <- Emod*(Amod+parv[6]*ABrad/2)*(1+(concA/parv[1])^parv[3])/(Aden*(parv[8]-parv[7]))
derivs[,9] <- Emod*(Bmod+parv[6]*ABrad/2)*(1+(concB/parv[2])^parv[4])/(Bden*(parv[9]-parv[7]))
derivs[Emodz,8:9] <- 0
derivs[,10] <- 1/Den-derivs[,8]*(parv[8]-parv[7])/(parv[10]-parv[7])-derivs[,9]*(parv[9]-parv[7])/(parv[10]-parv[7])
}
if (is.na(fixed[7])) { derivs[,7] <- 1-derivs[,8]-derivs[,9]-derivs[,10] }
}
return(cbind(Ei,derivs))
}
invertBRAIDrsm <- function(val,DA=NULL,DB=NULL,parv) {
if ((is.null(DA)&&is.null(DB))||(!is.null(DA)&&!is.null(DB))) {
stop("Exactly one of the inputs 'DA' and 'DB' must be NULL (the default).")
}
conc1 <- DA
conc2 <- DB
if (is.null(conc1)) {
if (length(val)!=length(conc2)) {
if (length(val)%%length(conc2)==0) { conc2 <- rep(conc2,times=length(val)/length(conc2)) }
else if (length(conc2)%%length(val)==0) { val <- rep(val,times=length(conc2)/length(val)) }
else { stop("The length of the longer of inputs 'val' and 'conc2' is not a multiple of the length of the shorter.") }
}
} else {
if (length(val)!=length(conc1)) {
if (length(val)%%length(conc1)==0) { conc1 <- rep(conc1,times=length(val)/length(conc1)) }
else if (length(conc1)%%length(val)==0) { val <- rep(val,times=length(conc1)/length(val)) }
else { stop("The length of the longer of inputs 'val' and 'conc1' is not a multiple of the length of the shorter.") }
}
}
concout <- rep(0,times=length(val))
incr <- sign(parv[10]-parv[7])
Eval <- (val-parv[7])/(parv[10]-val)
concout[which(incr*(val-parv[7])<=0)] <- 0
concout[which(incr*(val-parv[10])>=0)] <- Inf
cinds <- which(Eval>0)
Eval <- Eval[cinds]^(1/(parv[5]*sqrt(parv[3]*parv[4])))
if (is.null(conc1)) {
if (parv[9]==parv[10]) {
Bmod <- (conc2[cinds]/parv[2])^(sqrt(parv[4]/parv[3])/parv[5])
} else {
Bden <- 1+(1-((parv[9]-parv[7])/(parv[10]-parv[7])))*((conc2[cinds]/parv[2])^parv[4])
Bmod <- ((parv[9]-parv[7])/(parv[10]-parv[7]))*((conc2[cinds]/parv[2])^parv[4])
iinds <- which(is.infinite(Bden)|is.infinite(Bmod))
Bden[iinds] <- 1-(parv[9]-parv[7])/(parv[10]-parv[7])
Bmod[iinds] <- (parv[9]-parv[7])/(parv[10]-parv[7])
Bmod <- (Bmod/Bden)^(1/(parv[5]*sqrt(parv[3]*parv[4])))
}
if (parv[6]>=0) {
concout[cinds[which(Bmod>=Eval)]] <- 0
ccinds <- which(Bmod<Eval)
} else {
concout[cinds[which(Bmod>=Eval/(1-(parv[6]^2)/4))]] <- 0
ccinds <- which(Bmod<Eval/(1-(parv[6]^2)/4))
}
Bmod <- Bmod[ccinds]
Eval <- Eval[ccinds]
cinds <- cinds[ccinds]
Jval <- Eval+((parv[6]^2)/2-1)*Bmod
if (parv[6]!=0) { Jval <- Jval-parv[6]*sqrt(((parv[6]^2)/4-1)*Bmod^2+Eval*Bmod) }
Jval <- Jval^(parv[5]*sqrt(parv[3]*parv[4]))
if (parv[8]!=parv[10]) {
concout[cinds[which(Jval>=(parv[8]-parv[7])/(parv[10]-parv[8]))]] <- Inf
cinds <- cinds[which(Jval<(parv[8]-parv[7])/(parv[10]-parv[8]))]
Jval <- Jval[which(Jval<(parv[8]-parv[7])/(parv[10]-parv[8]))]
}
conc1c <- parv[1]*(Jval/((1+Jval)*(parv[8]-parv[7])/(parv[10]-parv[7])-Jval))^(1/parv[3])
concout[cinds] <- conc1c
} else {
if (parv[8]==parv[10]) {
Amod <- (conc1/parv[1])^(sqrt(parv[3]/parv[4])/parv[5])
} else {
Aden <- 1+(1-((parv[8]-parv[7])/(parv[10]-parv[7])))*((conc1/parv[1])^parv[3])
Amod <- ((parv[8]-parv[7])/(parv[10]-parv[7]))*((conc1/parv[1])^parv[3])
iinds <- which(is.infinite(Aden)|is.infinite(Amod))
Aden[iinds] <- 1-(parv[8]-parv[7])/(parv[10]-parv[7])
Amod[iinds] <- (parv[8]-parv[7])/(parv[10]-parv[7])
Amod <- (Amod/Aden)^(1/(parv[5]*sqrt(parv[3]*parv[4])))
}
if (parv[6]>=0) {
concout[cinds[which(Amod>=Eval)]] <- 0
ccinds <- which(Amod<Eval)
} else {
concout[cinds[which(Amod>=Eval/(1-(parv[6]^2)/4))]] <- 0
ccinds <- which(Amod<Eval/(1-(parv[6]^2)/4))
}
Amod <- Amod[ccinds]
Eval <- Eval[ccinds]
cinds <- cinds[ccinds]
Jval <- Eval+((parv[6]^2)/2-1)*Amod
if (parv[6]!=0) { Jval <- Jval-parv[6]*sqrt(((parv[6]^2)/4-1)*Amod^2+Eval*Amod) }
Jval <- Jval^(parv[5]*sqrt(parv[3]*parv[4]))
if (parv[9]!=parv[10]) {
concout[cinds[which(Jval>=(parv[9]-parv[7])/(parv[10]-parv[9]))]] <- Inf
cinds <- cinds[which(Jval<(parv[9]-parv[7])/(parv[10]-parv[9]))]
Jval <- Jval[which(Jval<(parv[9]-parv[7])/(parv[10]-parv[9]))]
}
conc2c <- parv[2]*(Jval/((1+Jval)*(parv[9]-parv[7])/(parv[10]-parv[7])-Jval))^(1/parv[4])
concout[cinds] <- conc2c
}
return(concout)
}
# (...,NA,NA,NA) : All are varying, EAB >= EA,EB
# (...,NA,NA,~) EAB=rmax(EA,EB) : EA and EB varying, EAB varies with maximum
# (...,NA,NA,~) EAB>rmax(EA,EB) : EA and EB varying below fixed EAB
# (...,NA,~,NA) EAB=EB>=EA : EA varying below EAB and EB varying together
# (...,NA,~,NA) EAB>EB;EAB>=EA : EA varying, EB fixed EAB varying above EA and EB
# (...,~,~,NA) EAB=EA=EB : All three varying together
# (...,~,~,NA) EAB=EA>EB : EAB and EA varying together above EB
# (...,~,~,NA) EAB>EA & EAB>EB : EAB varying above fixed EA and EB
# (...,NA,~,~) EAB=rmax(EA,EB) : EA varying, EB fixed, EAB varies with maximum
# (...,NA,~,~) EAB>EA;EAB>EB : EA varying below fixed EAB, EB fixed
# (...,~,~,~) : All three fixed
fitBRAIDrsm <- function(conc1,conc2,act,def=NULL,fixed=NULL,llims=NULL,ulims=NULL) {
# Fill in 'fixed' and 'def'
if (is.null(fixed)) { tfixed <- c(1,2,3,4,6,7,8,9) }
else { tfixed <- fixed }
if (length(tfixed)!=10 || length(which(is.na(tfixed)))==0) {
if (is.character(tfixed)) {
if (tfixed=="kappa1") { fixinds <- c(1,2,3,4,6,7,10) }
else if (tfixed=="kappa2") { fixinds <- c(1,2,3,4,6,7,8,9) }
else if (tfixed=="delta1") { fixinds <- c(1,2,3,4,5,7,10) }
else if (tfixed=="delta2") { fixinds <- c(1,2,3,4,5,7,8,9) }
else if (tfixed=="kappa3") { fixinds <- c(1,2,3,4,6,7,8,9,10) }
else if (tfixed=="delta3") { fixinds <- c(1,2,3,4,5,7,8,9,10) }
else if (tfixed=="ebraid") { fixinds <- 1:10 }
else { stop(sprintf("Unknown model name '%s'.",tfixed)) }
} else {
if (length(which(is.na(tfixed)))>0) { stop("'NA' values are only permitted in raw 'fixed' vectors.") }
fixinds <- sort(unique(round(tfixed)))
if (min(fixinds)<1 || max(fixinds)>10) { stop("Indexed 'fixed' vectors may only contain values from 1 to 10.") }
}
} else { fixinds <- which(is.na(tfixed)) }
npar <- length(fixinds)
if (npar==0) { stop("At least one parameter must be allowed to vary.") }
if (is.null(def)) {
inds <- which(conc1+conc2>0)
sdef <- dr_start(conc1[inds]+conc2[inds],act[inds])
odef <- c(exp(sdef[4]),exp(sdef[4]),sdef[1],sdef[1],1,0,sdef[2],sdef[3],sdef[3],sdef[3])
if (length(tfixed)!=10 || length(which(is.na(tfixed)))==0) {
tfixed <- odef
tfixed[fixinds] <- NA
} else {
odef[which(!is.na(tfixed))] <- tfixed[which(!is.na(tfixed))]
if (odef[10]>odef[7] && !is.na(tfixed[10])) { odef[8:10] <- pmin(odef[8:10],tfixed[10]) }
else if (!is.na(tfixed[10])) { odef[8:10] <- pmax(odef[8:10],tfixed[10]) }
}
} else {
if (length(tfixed)!=10 || length(which(is.na(tfixed)))==0) {
if (length(def)==10) {
odef <- def
tfixed <- odef
tfixed[fixinds] <- NA
} else { stop("When using index vectors or model names for 'fixed', parameter 'def' must specify all 10 values.") }
} else {
if (length(def)==10) {
odef <- def
tfixed <- odef
tfixed[fixinds] <- NA
} else if (length(def)==npar) {
odef <- tfixed
odef[fixinds] <- def
} else { stop("Input 'def' must have as many values as free parameters or all parameter values (length 10).") }
}
}
fixed <- tfixed
isIncr <- sign(odef[10]-odef[7])
# medAct <- isIncr*min(max(isIncr*median(act),isIncr*odef[7]),min(isIncr*odef[8:10]))
medAct <- isIncr*(isIncr*odef[7]+min(isIncr*odef[8:10]))/2
# Fill in lower limit values
if (is.null(llims)) { llims <- rep(NA,times=10) }
else {
if (length(llims)<10 && length(llims)!=npar) { stop("Improper length for lower limit vector.") }
else if (length(llims)<10) {
tllims <- rep(NA,times=10)
tllims[which(is.na(fixed))] <- llims
llims <- tllims
}
}
full_llims <- c(min(exp(-12),odef[1]/exp(3)),min(exp(-12),odef[2]/exp(3)),1/10,1/10,1/10,-1.999,0,0,0,0)
if (isIncr>0) { full_llims[7:8] <- c(2*odef[7]-odef[10],medAct+0.001) }
else if (isIncr<0) { full_llims[7:8] <- c(medAct+0.001,2*odef[10]-odef[7]) }
else { stop("Default limits must cover a non-zero width range.") }
full_llims[9:10] <- full_llims[8]
full_llims[which(!is.na(llims))] <- pmax(full_llims[which(!is.na(llims))],llims[which(!is.na(llims))])
full_llims[1] <- max(full_llims[1],exp(2*min(log(conc1[conc1>0]))-max(log(conc1[conc1>0]))))
full_llims[2] <- max(full_llims[2],exp(2*min(log(conc2[conc2>0]))-max(log(conc2[conc2>0]))))
odef[which(is.na(fixed))] <- pmax(odef,full_llims)[which(is.na(fixed))]
# Fill in upper limit values
if (is.null(ulims)) { ulims <- rep(NA,times=10) }
else {
if (length(ulims)<10 && length(ulims)!=npar) { stop("Improper length for upper limit vector.") }
else if (length(ulims)<10) {
tulims <- rep(NA,times=10)
tulims[which(is.na(fixed))] <- ulims
ulims <- tulims
}
}
full_ulims <- c(max(1,odef[1]*exp(3)),max(1,odef[2]*exp(3)),10,10,10,100,0,0,0,0)
if (isIncr>0) { full_ulims[7:8] <- c(medAct-0.001,2*odef[10]-odef[7]) }
else { full_ulims[7:8] <- c(2*odef[7]-odef[10],medAct-0.001) }
full_ulims[9:10] <- full_ulims[8]
full_ulims[which(!is.na(ulims))] <- pmin(full_ulims[which(!is.na(ulims))],ulims[which(!is.na(ulims))])
full_ulims[1] <- min(full_ulims[1],exp(2*max(log(conc1[conc1>0]))-min(log(conc1[conc1>0]))))
full_ulims[2] <- min(full_ulims[2],exp(2*max(log(conc2[conc2>0]))-min(log(conc2[conc2>0]))))
odef[which(is.na(fixed))] <- pmin(odef,full_ulims)[which(is.na(fixed))]
# Prepare parameter values and bounds for optimization
tdef <- odef
tdef[1:5] <- log(tdef[1:5])
full_llims[1:5] <- log(full_llims[1:5])
full_ulims[1:5] <- log(full_ulims[1:5])
if (is.na(fixed[10])) {
if (is.na(fixed[8])||is.na(fixed[9])) {
if (!is.na(fixed[8]) && odef[8]==odef[10]) {
tdef[10] <- tdef[10]-tdef[9]
if (isIncr>0) {
full_llims[10] <- 0
full_ulims[10] <- full_ulims[10]-tdef[9]
} else {
full_llims[10] <- full_llims[10]-tdef[7]
full_ulims[10] <- 0
}
} else if (!is.na(fixed[9]) && odef[9]==odef[10]) {
tdef[10] <- tdef[10]-tdef[8]
if (isIncr>0) {
full_llims[10] <- 0
full_ulims[10] <- full_ulims[10]-tdef[8]
} else {
full_llims[10] <- full_llims[10]-tdef[8]
full_ulims[10] <- 0
}
} else {
tdef[10] <- tdef[10]-isIncr*max(isIncr*tdef[8:9])
if (isIncr>0) {
full_llims[10] <- 0
full_ulims[10] <- full_ulims[10]-isIncr*max(isIncr*tdef[8:9])
} else {
full_llims[10] <- full_llims[10]-isIncr*max(isIncr*tdef[8:9])
full_ulims[10] <- 0
}
}
} else {
if (odef[8]!=odef[10]) {
if (isIncr>0) { full_llims[10] <- max(full_llims[10],odef[8]) }
else { full_ulims[10] <- min(full_ulims[10],odef[8]) }
}
if (odef[9]!=odef[10]) {
if (isIncr>0) { full_llims[10] <- max(full_llims[10],odef[9]) }
else { full_ulims[10] <- min(full_ulims[10],odef[9]) }
}
}
} else if (isIncr*odef[10]>max(isIncr*odef[8:9])) {
if (isIncr>0) { full_ulims[8:9] <- pmin(full_ulims[8:9],odef[10]) }
else { full_llims[8:9] <- pmax(full_llims[8:9],odef[10]) }
}
llims <- full_llims[which(is.na(fixed))]
ulims <- full_ulims[which(is.na(fixed))]
parv2fullpar <- function(parv) {
fullpar <- tdef
fullpar[which(is.na(fixed))] <- parv
fullpar[1:5] <- exp(fullpar[1:5])
if (is.na(fixed[10])) {
if (is.na(fixed[8])||is.na(fixed[9])) {
if (!is.na(fixed[8]) && odef[8]==odef[10]) {
fullpar[10] <- fullpar[10]+fullpar[9]
fullpar[8] <- fullpar[10]
} else if (!is.na(fixed[9]) && odef[9]==odef[10]) {
fullpar[10] <- fullpar[10]+fullpar[8]
fullpar[9] <- fullpar[10]
} else { fullpar[10] <- fullpar[10]+isIncr*max(isIncr*fullpar[8:9]) }
} else {
if (odef[8]==odef[10]) { fullpar[8] <- fullpar[10] }
if (odef[9]==odef[10]) { fullpar[9] <- fullpar[10] }
}
} else if (is.na(fixed[8]) && is.na(fixed[9])) {
if (odef[10]==isIncr*max(isIncr*odef[8:9])) { fullpar[10] <- isIncr*max(isIncr*fullpar[8:9]) }
} else if (is.na(fixed[8]) && ((odef[8]==odef[10]&&isIncr*fullpar[8]>=isIncr*fullpar[9])
|| isIncr*fullpar[8]>isIncr*odef[10])) { fullpar[10] <- fullpar[8] }
else if (is.na(fixed[9]) && ((odef[9]==odef[10]&&isIncr*fullpar[9]>=isIncr*fullpar[8])
|| isIncr*fullpar[9]>isIncr*odef[10])) { fullpar[10] <- fullpar[9] }
return(fullpar)
}
parvRSMErr <- function(parv) {
fullpar <- parv2fullpar(parv)
ei <- evalBRAIDrsm(conc1,conc2,fullpar)
err <- ei-act
return(sum(err^2))
}
parvRSMDerivs <- function(parv) {
fullpar <- parv2fullpar(parv)
ederivs <- evalBRAIDrsm_int(conc1,conc2,fullpar,fixed=fixed,calcderivs=TRUE)
err <- ederivs[,1]-act
for (i in 1:5) { ederivs[,i+1] <- fullpar[i]*ederivs[,i+1] }
if (is.na(fixed[10]) && (is.na(fixed[8])||is.na(fixed[9]))) {
if (!is.na(fixed[8]) && odef[8]==odef[10]) { ederivs[,11] <- ederivs[,11]-ederivs[,10] }
else if (!is.na(fixed[9]) && odef[9]==odef[10]) { ederivs[,11] <- ederivs[,11]-ederivs[,9] }
else {
if (isIncr*fullpar[8]>=isIncr*fullpar[9]) { ederivs[,11] <- ederivs[,11]-ederivs[,9] }
else { ederivs[,11] <- ederivs[,11]-ederivs[,10] }
}
}
oderivs <- rep(0,times=10)
for (i in 1:10) { oderivs[i] <- 2*sum(err*ederivs[,i+1]) }
oderivs <- oderivs[which(is.na(fixed))]
return(oderivs)
}
escale <- max(abs(odef[5]-odef[9]),abs(odef[6]-odef[9]),abs(odef[10]-odef[9]))
parscale <- c(1,1,1,1,1,1,0.1*escale,0.1*escale,0.1*escale,0.1*escale)
control <- list(maxit=500,parscale=parscale[which(is.na(fixed))])
nls <- optim(tdef[which(is.na(fixed))],parvRSMErr,parvRSMDerivs,method="L-BFGS-B",
lower=llims,upper=ulims,hessian=FALSE,control=control)
nls$fullpar <- parv2fullpar(nls$par)
nls$fixed <- fixed
nls$odef <- odef
full_llims <- parv2fullpar(llims)
llims <- full_llims[which(is.na(fixed))]
full_ulims <- parv2fullpar(ulims)
ulims <- full_ulims[which(is.na(fixed))]
nls$mlims <- rbind(llims,ulims)
return(nls)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.