Description Usage Arguments Details Value Author(s) See Also Examples
Implements drc nonlinear functions into the nlme framework for mixed effects dose-response modeling.
1 2 |
form |
Formula describing the dose-response relationship |
curveid |
Formula with parameter names on the left hand side (divided by +) and a column name in data, denoting a factor, to estimate separate parameters per factor-level. If NULL only fixed effects for a single curve will be estimated. |
data |
a data.frame object |
fct |
a function of the drc package |
random |
a list or one-sided formula describing the random effects |
correlation |
additional corClasses object |
weights |
additional varClasses object |
control |
list with nlme control arguments |
start |
optional list with initial values for the fixed components. If NULL the initial values will be found automatically. |
REML |
logical value. If TRUE the model is fit by maximizing the restricted log-likelihood, if FALSE log-likelihood is maximized. |
An application of medrm is shown on the help pages of data broccoli
. EDx and selectivity indices can be calculated with functions ED
and EDcomp
. Model-averaged ED can be computed by function mmaED
.
Marginal EDx and corresponding function are implemented in functions EDmarg
, using numerical integration conditional on the etimated variance components.
An object of class medrc
Daniel Gerhard
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 | library(nlme)
data(DNase)
DNase$Run <- factor(DNase$Run, levels=1:11)
ggplot(DNase, aes(y=density, x=conc, colour=Run)) + geom_point() + coord_trans(x="log")
#############
# fit a 5-parameter log-logistic model
# with random asymptotes, steepness and inflection point location
## starting values need to be provided:
# estimate a set of parameters for each Run
mf <- glsdrm(density ~ conc, curveid=b + c + d + e + f ~ Run,
fct=LL.5(), data=DNase)
# curve specific difference to the parameter average
cmat <- matrix(coefficients(mf), ncol=5)
m <- apply(cmat, 2, mean) # fixed effect starting values
names(m) <- letters[2:6]
rmat <- t(apply(cmat, 1, function(x) x-m)) # random effect starting values
rownames(rmat) <- levels(DNase$Run)
colnames(rmat) <- letters[2:6]
# mixed model fit
mm <- medrm(density ~ conc, fct=LL.5(), data=DNase,
random=b + c + d + e ~ 1 | Run,
start=list(fixed=m[c("b","c","d","e","f")], Run=rmat))
# medrc methods:
plot(mm, logx=TRUE, ndose=100, ranef=TRUE)
print(mm)
VarCorr(mm)
ranef(mm)
coef(mm)
vcov(mm)
### effective dose estimation ED25, ED50, ED75
# based on the fixed effect estimates
# conditional on random effects of 0
ED(mm, c(25, 50, 75), interval="tfls")
# predictions for Run 2
ndata <- data.frame(conc=exp(seq(log(0.1),log(13), length=10)), Run="2")
predict(mm, newdata=ndata)
##############################
# fitting a reduced model
# without random slope ()
mmr <- medrm(density ~ conc, fct=LL.5(), data=DNase,
random=c + d + e ~ 1 | Run,
start=list(fixed=m[c("b","c","d","e","f")], Run=rmat))
# estimating marginal ED25, ED50, ED75
# (by numerical integration conditional on estimated variance-components)
EDmarg(mmr, c(25, 50, 75), interval="tfls")
# comparing to Run-specific ED
ED(mmr, c(25, 50, 75), interval="tfls")
# marginal prediction
ndata <- data.frame(conc=exp(seq(log(0.1),log(13), length=10)))
predict(mmr, type="marginal", newdata=ndata)
#############
#############
# fit a 4-parameter Weibull model
# with random asymptotes and inflection point location
## starting values:
mf <- glsdrm(density ~ conc, curveid=b + c + d + e ~ Run,
fct=W1.4(), data=DNase)
cmat <- matrix(coefficients(mf), ncol=4)
m <- apply(cmat, 2, mean) # fixed effect starting values
names(m) <- letters[2:5]
rmat <- t(apply(cmat, 1, function(x) x-m)) # random effect starting values
rownames(rmat) <- levels(DNase$Run)
colnames(rmat) <- letters[2:5]
# mixed model fit
mmw <- medrm(density ~ conc, fct=W1.4(), data=DNase,
random=c + d + e ~ 1 | Run,
start=list(fixed=m[c("b","c","d","e")], Run=rmat))
plot(mmw, logx=TRUE, ndose=100, ranef=TRUE)
# AIC comparison of LL.5 and W1.4 model
AIC(mm, mmw)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.