

   MMuullttiivvaarriiaattee eelllliippttiiccaallllyy--ccoonnttoouurreedd rreeppeeaatteedd mmeeaassuurreemmeennttss
   mmooddeellss ffoorr lliinneeaarr aanndd nnoonnlliinneeaarr cchhaannggeess oovveerr ttiimmee iinn tthhee
   pprreesseennccee ooff ttiimmee--vvaarryyiinngg ccoovvaarriiaatteess aanndd wwiitthh AARR((11)) aanndd ttwwoo
   lleevveellss ooff vvaarriiaannccee ccoommppoonneennttss

        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)

   AArrgguummeennttss::

   response: A list of two or three column matrices with
             response values, times, and possibly nesting cate-
             gories, 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 polyno-
             mial in time is automatically fitted; (2) `logis-
             tic' for a four-parameter logistic growth curve;
             (3) `pkpd' for a first-order one-compartment phar-
             macokinetic model.  Otherwise, set this to a func-
             tion 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 can-
             not 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 neces-
             sary 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 base-
             line 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 vari-
             ables must be binary (or factor variables) as dif-
             ferent 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 func-
             tion.  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 can-
             not 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 increas-
             ing integers. This option should always be speci-
             fied 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 interac-
             tions between them and the polynomial in time in
             the `linear model'.

   transform: Transformation of the response variable: `iden-
             tity', `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: `expo-
             nential' 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 other-
             wise; `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 func-
             tion 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 con-
             tains two initial values, an additive constant is
             included: pvar(1)+pvar(2)*v(mu). Otherwise, either
             a function or a formula beginning with ~, specify-
             ing either a linear regression function in the
             Wilkinson and Rogers notation or a general func-
             tion 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 vari-
             ance components, depending on the number of levels
             of nesting.

        par: If supplied, an initial estimate for the autocor-
             relation 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 dec-
             imals, 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'.

   DDeessccrriippttiioonn::

        `elliptic' fits a special case of the multivariate
        elliptically-contoured distribution, called the multi-
        variate 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 ran-
        dom 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 `fin-
        terp'.)

   VVaalluuee::

        A list of class `elliptic' is returned.

   AAuutthhoorr((ss))::

        J.K. Lindsey

   SSeeee AAllssoo::

        `carma', `finterp', `gar', `gettvc', `glmm', `gnlmm',
        `gnlr', `kalseries', `potthoff', `read.list',
        `restovec', `rmna', `tcctomat', `tvctomat'.

   EExxaammpplleess::

        # 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)

