elliptic(response, model="linear", times=NULL, dose=NULL, ccov=NULL, tvcov=NULL, nest=NULL, torder=0, interaction=NULL, transform="identity", link="identity", autocorr="exponential", pell=NULL, preg=rep(1,4), pvar=var(y), varfn=NULL, pre=NULL, par=NULL, delta=NULL, envir=sys.frame(sys.parent()), print.level=0, gradtol=0.00001, typsiz=abs(theta), stepmax=10*sqrt(theta%*%theta), steptol=0.00001, iterlim=100, ndigit=10, fscale=1)
response
|
A list of two or three column matrices with response
values, times, and possibly nesting categories, for each individual,
one matrix or dataframe of response values, or an object of class,
response (created by restovec) or repeated (created by
rmna).
|
model
|
The model to be fitted for the location. Builtin choices
are (1) linear for linear models with time-varying covariate; if
torder > 0, a polynomial in time is automatically fitted; (2)
logistic for a four-parameter logistic growth curve; (3)
pkpd for a first-order one-compartment pharmacokinetic model.
Otherwise, set this to a function of the parameters or a formula
beginning with ~, specifying either a linear regression function for
the location parameter in the Wilkinson and Rogers notation or a
general function with named unknown parameters that describes the
location, returning a vector the same length as the number of
observations, in which case ccov and tvcov cannot be used.
|
times
|
When response is a matrix, a vector of possibly
unequally spaced times when they are the same for all individuals or a
matrix of times. Not necessary if equally spaced. Ignored if response
has class, response or repeated.
|
dose
|
A vector of dose levels for the pkpd model, one per
individual.
|
ccov
|
A vector or matrix containing time-constant baseline
covariates with one line per individual, a model formula using
vectors of the same size, or an object of class, tccov (created by
tcctomat). For the pkpd and logistic
models, all variables must be binary (or factor variables) as
different values of all parameters are calculated for all combinations
of these variables (except for the logistic model when a time-varying
covariate is present). Cannot be used when model is a function.
Ignored if response has class, repeated.
|
tvcov
|
A list of vectors or matrices with time-varying
covariates for each individual (one column per variable), a matrix or
dataframe of such covariate values (if only one covariate), or an
object of class, tvcov (created by tvctomat). If times
are not the same as for responses, the list can be created with
gettvc. Only one time-varying covariate is allowed
except for the linear model; if more are required, set
model equal to the appropriate mean function. This argument
cannot be used when model is a function and is ignored if response has
class, repeated.
|
nest
|
When response is a matrix, a vector of length equal
to the number of responses per individual indicating which responses
belong to which nesting category. Categoriess must be consecutive
increasing integers. This option should always be specified if nesting
is present. Ignored if response has class, repeated.
|
torder
|
When the linear model is chosen, order of the
polynomial in time to be fitted.
|
interaction
|
Vector of length equal to the number of
time-constant covariates, giving the levels of interactions between
them and the polynomial in time in the linear model.
|
transform
|
Transformation of the response variable: identity,
exp, square, sqrt, or log.
|
link
|
Link function for the location: identity, exp,
square, sqrt, or log. For the linear model,
if not the identity, initial estimates of the regression
parameters must be supplied (intercept, polynomial in time,
time-constant covariates, time-varying covariates, in that order).
|
autocorr
|
The form of the autocorrelation function:
exponential is the usual rho^|t_i-t_j|; gaussian is
rho^((t_i-t_j)^2); cauchy is 1/(1+rho(t_i-t_j)^2);
spherical is ((|t_i-t_j|rho)^3-3|t_i-t_j|rho+2)/2 for
|t_i-t_j|<=1/rho and zero otherwise; IOU is the integrated
Ornstein-Uhlenbeck process, (2rho min(t_i,t_j)+exp(-rho t_i)
+exp(-rho t_j)-1 -exp(rho|ti-t_j|))/2rho^3.
|
pell
| Initial estimate of the power parameter of the multivariate elliptically-contoured distribution. If missing, the multivariate normal distribution is used. |
preg
|
Initial parameter estimates for the regression model.
Only required for linear model if the link is not the
identity or a variance function is fitted.
|
pvar
|
Initial parameter estimate for the variance. If more than
one value is provided, the log variance depends on a polynomial in
time. With the pkpd model, if four values are supplied, a
nonlinear regression for the variance is fitted.
|
varfn
|
The builtin variance function has the variance
proportional to a function of the location: pvar*v(mu) = identity
or square. If pvar contains two initial values, an additive
constant is included: pvar(1)+pvar(2)*v(mu). Otherwise, either a
function or a formula beginning with ~, specifying either a linear
regression function in the Wilkinson and Rogers notation or a general
function with named unknown parameters for the log variance can be
supplied, yielding a vector the same length as the number of
observations.
|
pre
| Zero, one or two parameter estimates for the variance components, depending on the number of levels of nesting. |
par
| If supplied, an initial estimate for the autocorrelation parameter. |
delta
| Scalar or vector giving the unit of measurement for each response value, set to unity by default. For example, if a response is measured to two decimals, delta=0.01. Ignored if response has class, response or repeated. |
envir
| Environment in which model formulae are to be interpreted or a data object of class, repeated, tccov, or tvcov. |
others
|
Arguments controlling nlm.
|
elliptic fits a special case of the multivariate
elliptically-contoured distribution, called the multivariate power
exponential distribution. It includes the multivariate normal
(power=1), the multivariate Laplace (power=0.5), and the multivariate
uniform (power -> infinity) distributions as special cases.
For clustered (non-longitudinal) data, where only random effects will be fitted, the `times' may be any strictly increasing sequence distinguishing the responses on an individual.
It is designed to fit linear and nonlinear models with time-varying covariates observed at arbitrary time points. A continuous-time AR(1) and zero, one, or two levels of nesting can be handled.
Nonlinear regression models can be supplied as formulae where
parameters are unknowns. Factor variables cannot be used and
parameters must be scalars. (See finterp.)
elliptic is returned.carma, finterp, gar,
gettvc, glmm, gnlmm,
gnlr, kalseries, potthoff,
read.list, restovec, rmna,
tcctomat, tvctomat.
# linear models
y <- matrix(rnorm(40),ncol=5)
x1 <- gl(2,4)
x2 <- gl(2,1,8)
# independence with time trend
elliptic(y, ccov=~x1, torder=2)
# AR(1)
elliptic(y, ccov=~x1, torder=2, par=0.1)
elliptic(y, ccov=~x1, torder=3, interact=3, par=0.1)
# random intercept
elliptic(y, ccov=~x1+x2, interact=c(2,0), torder=3, pre=2)
#
# nonlinear models
times <- rep(1:20,2)
dose <- c(rep(2,20),rep(5,20))
mu <- function(p) exp(p[1]-p[3])*(dose/(exp(p[1])-exp(p[2]))*
(exp(-exp(p[2])*times)-exp(-exp(p[1])*times)))
shape <- function(p) exp(p[1]-p[2])*times*dose*exp(-exp(p[1])*times)
conc <- matrix(rnorm(40,mu(log(c(1,0.3,0.2))),sqrt(shape(log(c(0.1,0.4))))),
ncol=20,byrow=T)
conc[,2:20] <- conc[,2:20]+0.5*(conc[,1:19]-matrix(mu(log(c(1,0.3,0.2))),
ncol=20,byrow=T)[,1:19])
conc <- ifelse(conc>0,conc,0.01)
# with builtin function
# independence
elliptic(conc, model="pkpd", preg=log(c(0.5,0.4,0.1)), dose=c(2,5))
# AR(1)
elliptic(conc, model="pkpd", preg=log(c(0.5,0.4,0.1)), dose=c(2,5),
par=0.1)
# add variance function
elliptic(conc, model="pkpd", preg=log(c(0.5,0.4,0.1)), dose=c(2,5),
par=0.1, varfn=shape, pvar=log(c(0.5,0.2)))
# multivariate elliptical distribution
elliptic(conc, model="pkpd", preg=log(c(0.5,0.4,0.1)), dose=c(2,5),
par=0.1, varfn=shape, pvar=log(c(0.5,0.2)), pell=1)
# or equivalently with user-specified function
# independence
elliptic(conc, model=mu, preg=log(c(0.5,0.4,0.1)))
# AR(1)
elliptic(conc, model=mu, preg=log(c(0.5,0.4,0.1)), par=0.1)
# add variance function
elliptic(conc, model=mu, preg=log(c(0.5,0.4,0.1)), par=0.1,
varfn=shape, pvar=log(c(0.5,0.2)))
# multivariate elliptical distribution
elliptic(conc, model=mu, preg=log(c(0.5,0.4,0.1)), par=0.1,
varfn=shape, pvar=log(c(0.5,0.2)), pell=1)
# or with user-specified formula
# independence
elliptic(conc, model=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
preg=list(absorption=log(0.5),elimination=log(0.4),
volume=log(0.1)))
# AR(1)
elliptic(conc, model=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
preg=list(absorption=log(0.5),elimination=log(0.4),volume=log(0.1)),
par=0.1)
# add variance function
elliptic(conc, model=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
preg=list(absorption=log(0.5),elimination=log(0.4),volume=log(0.1)),
varfn=~exp(b1-b2)*times*dose*exp(-exp(b1)*times),
par=0.1, pvar=list(b1=log(0.5),b2=log(0.2)))
# multivariate elliptical distribution
elliptic(conc, model=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
preg=list(absorption=log(0.5),elimination=log(0.4),volume=log(0.1)),
varfn=~exp(b1-b2)*times*dose*exp(-exp(b1)*times),
par=0.1, pvar=list(b1=log(0.5),b2=log(0.2)), pell=1)
#
# generalized logistic regression with square-root transformation
# and square link
times <- rep(seq(10,200,by=10),2)
mu <- function(p) {
yinf <- exp(p[2])
yinf*(1+((yinf/exp(p[1]))^p[4]-1)*exp(-yinf^p[4]
*exp(p[3])*times))^(-1/p[4])}
y <- matrix(rnorm(40,sqrt(mu(c(2,1.5,0.05,-2))),0.05)^2,ncol=20,byrow=T)
y[,2:20] <- y[,2:20]+0.5*(y[,1:19]-matrix(mu(c(2,1.5,0.05,-2)),
ncol=20,byrow=T)[,1:19])
y <- ifelse(y>0,y,0.01)
# with builtin function
# independence
elliptic(y, model="logistic", preg=c(2,1,0.1,-1), trans="sqrt",
link="square")
# the same model with AR(1)
elliptic(y, model="logistic", preg=c(2,1,0.1,-1), trans="sqrt",
link="square", par=0.4)
# the same model with AR(1) and one component of variance
elliptic(y, model="logistic", preg=c(2,1,0.1,-1),
trans="sqrt", link="square", pre=1, par=0.4)
# or equivalently with user-specified function
# independence
elliptic(y, model=mu, preg=c(2,1,0.1,-1), trans="sqrt",
link="square")
# the same model with AR(1)
elliptic(y, model=mu, preg=c(2,1,0.1,-1), trans="sqrt",
link="square", par=0.4)
# the same model with AR(1) and one component of variance
elliptic(y, model=mu, preg=c(2,1,0.1,-1),
trans="sqrt", link="square", pre=1, par=0.4)
# or equivalently with user-specified formula
# independence
elliptic(y, model=~exp(yinf)*(1+((exp(yinf-y0))^b4-1)*
exp(-exp(yinf*b4+b3)*times))^(-1/b4),
preg=list(y0=2,yinf=1,b3=0.1,b4=-1), trans="sqrt", link="square")
# the same model with AR(1)
elliptic(y, model=~exp(yinf)*(1+((exp(yinf-y0))^b4-1)*
exp(-exp(yinf*b4+b3)*times))^(-1/b4),
preg=list(y0=2,yinf=1,b3=0.1,b4=-1), trans="sqrt",
link="square", par=0.1)
# add one component of variance
elliptic(y, model=~exp(yinf)*(1+((exp(yinf-y0))^b4-1)*
exp(-exp(yinf*b4+b3)*times))^(-1/b4),
preg=list(y0=2,yinf=1,b3=0.1,b4=-1),
trans="sqrt", link="square", pre=1, par=0.1)
#
# multivariate elliptical distribution for outliers
y <- matrix(rcauchy(40,mu(c(2,1.5,0.05,-2)),0.05),ncol=20,byrow=T)
y[,2:20] <- y[,2:20]+0.5*(y[,1:19]-matrix(mu(c(2,1.5,0.05,-2)),
ncol=20,byrow=T)[,1:19])
y <- ifelse(y>0,y,0.01)
# first with normal distribution
elliptic(y, model="logistic", preg=c(1,1,0.1,-1))
elliptic(y, model="logistic", preg=c(1,1,0.1,-1), par=0.5)
# then elliptic
elliptic(y, model="logistic", preg=c(1,1,0.1,-1), pell=1)
elliptic(y, model="logistic", preg=c(1,1,0.1,-1), par=0.5, pell=1)