# R/mixture.R In drc: Analysis of Dose-Response Curves

```"mixture" <- function(object, model = c("CA", "Hewlett", "Voelund"), start, startm, control = drmc())
{
model <- match.arg(model)

if (!is.null(object\$"call"\$"pmodels"[[2]]))
{
bListEl <- as.formula(object\$"call"\$"pmodels"[[2]])
} else {
stop("'collapse' argument should be a formula")
}

## Constructing new collapse argument ... in common for Hewlett and Voelund
curveid <- as.character(object\$"call"\$"curve")
eListEl1 <- as.formula(paste("~ I(1/(", curveid, "/100))-1", sep = ""))
eListEl2 <- as.formula(paste("~ I(1/(1-", curveid, "/100))-1", sep = ""))

blstr <- paste(bListEl, collapse = "")
el1str <- paste(eListEl1, collapse = "")
el2str <- paste(eListEl2, collapse = "")
#
#    print(bListEl)
#    print(eListEl1)
#    print(eListEl2)

fct <- object\$"fct"
if (model == "CA")
{
if (length(fct\$names) == 2)  # LL.2
{
mixtfct <- hewlett(fixed = c( NA, 0, 1, NA, NA, 1))
collapseNew2 <- list(bListEl, eListEl1, eListEl2)
pmText <- paste(blstr, el1str, el2str, sep = ", ")

eName <- fct\$names[2]
noLim <- 0  # number of c and d parameters
} else if (length(fct\$names) == 3)  # l3 function
{
mixtfct <- hewlett(fixed = c( NA, 0, NA, NA, NA, 1))
collapseNew2 <- list(bListEl, ~1, eListEl1, eListEl2)
pmText <- paste(blstr, "~1", el1str, el2str, sep = ", ")

eName <- fct\$names[3]
noLim <- 1  # number of c and d parameters
} else if (length(fct\$names) == 4)  # l4 function
{
mixtfct <- hewlett(fixed = c( NA, NA, NA, NA, NA, 1))
collapseNew2 <- list(bListEl, ~1, ~1, eListEl1, eListEl2)
pmText <- paste(blstr, "~1", "~1", el1str, el2str, sep = ", ")

eName <- fct\$names[4]
noLim <- 2  # number of c and d parameters
} else {
stop("Does not work for LL.5()")
}

if (missing(startm)) {startm <- NULL}

class(mixtfct) <- "CA"  # otherwise it would be "Hewlett"
mixtfct\$"name" <- "ca"  # overruling the default name "hewlett"
}

if (model == "Hewlett")
{
if (length(fct\$names) == 2)  # LL.2
{
mixtfct <- hewlett(fixed = c(NA, 0, 1, NA, NA, NA))
collapseNew2 <- list(bListEl, eListEl1, eListEl2, ~1)
pmText <- paste(blstr, el1str, el2str, "~1", sep = ", ")

eName <- fct\$names[2]
noLim <- 0  # number of c and d parameters
} else if (length(fct\$names) == 3)  # LL.3
{
mixtfct <- hewlett(fixed = c(NA, 0, NA, NA, NA, NA))
collapseNew2 <- list(bListEl, ~1, eListEl1, eListEl2, ~1)
pmText <- paste(blstr, "~1", el1str, el2str, "~1", sep = ", ")

eName <- fct\$names[3]
noLim <- 1  # number of c and d parameters
} else if (length(fct\$names) == 4)  # LL.4
{
mixtfct <- hewlett()
collapseNew2 <- list(bListEl, ~1, ~1, eListEl1, eListEl2, ~1)
pmText <- paste(blstr, "~1", "~1", el1str, el2str, "~1", sep = ", ")

eName <- fct\$names[4]
noLim <- 2  # number of c and d parameters
} else {
stop("Does not work for LL.5()")
}
if (missing(startm)) {startm <- 1}
}

if (model == "Voelund")
{
if (length(fct\$names) == 2)  # LL.2
{
mixtfct <- voelund(fixed = c( NA, 0, 1, NA, NA, NA, NA))
collapseNew2 <- list(bListEl, eListEl1, eListEl2, ~1, ~1)
pmText <- paste(blstr, el1str, el2str, "~1", "~1", sep = ", ")

eName <- fct\$names[2]
noLim <- 0  # number of c and d parameters
} else if (length(fct\$names) == 3)  # LL.3
{
mixtfct <- voelund(fixed = c( NA, 0, NA, NA, NA, NA, NA))
collapseNew2 <- list(bListEl, ~1, eListEl1, eListEl2, ~1, ~1)
pmText <- paste(blstr, "~1", el1str, el2str, "~1", "~1", sep = ", ")

eName <- fct\$names[3]
noLim <- 1  # number of c and d parameters
} else if (length(fct\$names) == 4)  # LL.4
{
mixtfct <- voelund()
collapseNew2 <- list(bListEl, ~1, ~1, eListEl1, eListEl2, ~1, ~1)
pmText <- paste(blstr, "~1", "~1", el1str, el2str, "~1", "~1", sep = ", ")

eName <- fct\$names[4]
noLim <- 2  # number of c and d parameters
} else {
stop("Does not work for LL.5()")
}
if (missing(startm)) {startm <- c(3, 0.3)}
}

#    assign("collapseNew2", collapseNew2, envir = .GlobalEnv)

## Checking if levels 0 and 100 are present
assayNo <- object\$"dataList"\$"curveid"
if (all(regexpr("0", as.character(unique(assayNo))) < 0 ))
{
stop("Level 0 is missing")
}
if (all(regexpr("100", as.character(unique(assayNo))) < 0 ))
{
stop("Level 100 is missing")
}

## Constructing starting values
sv <- coef(object)

parNames1 <- object\$"parNames"[[1]]
parNames2 <- object\$"parNames"[[3]]
#    eNames <- as.character(parNames2[regexpr(paste(eName, ":", sep=""), parNames1, fixed = TRUE) > 0])
#    eInd <- grep(paste(eName, ":", sep = ""), parNames1)
#    eNames <- as.character(parNames2[eInd])

#    pos0 <- match(paste("factor(", curveid, ")0", sep = ""), eNames)
#    pos1 <- match(paste("factor(", curveid, ")100", sep = ""), eNames)
#    if (is.na(pos0)) {pos0 <- 1}  # it is the intercept
#    if (is.na(pos1)) {pos1 <- 1}
#
#    noED50 <- length(eNames)
#    noB <- length(coef(object)) - noED50 - noLim
#    sv2 <- sv[c(1:(noB + noLim), noB + noLim + pos0, noB + noLim + pos1)]
#    sv2[noB+noLim+2] <- sv2[noB+noLim+1] + sv2[noB+noLim+2]

if (missing(start))
{
sv2 <- c(sv[-grep(paste(eName, ":", sep = ""), parNames1)],
sv[grep(paste(eName, ":", "100", sep = ""), parNames1)],
sv[grep(paste(eName, ":", "0", sep = ""), parNames1)])
sv3 <- as.vector(c(sv2, startm))  # removing the element names
} else {
sv3 <- start
}

#    print(mixtfct)
#    print(collapseNew2)
#    print(sv3)
mModel <- update(object, fct = mixtfct, pmodels = collapseNew2, start = sv3, control = control)

#    rm(collapseNew2, envir = .GlobalEnv)

mModel\$deviance <- object\$"fit"\$"value"

mModel\$"anova"\$"test" <- "F"
mModel\$"anova"\$"anovaFit" <- object\$"anova"\$"anovaFit"  # model1
mModel\$"text" <- paste(model, "model")
mModel\$"pmodelsText" <- pmText

return(mModel)
}
```

## Try the drc package in your browser

Any scripts or data that you put into this service are public.

drc documentation built on May 29, 2017, 11:31 a.m.