

   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 aanndd SSttuuddeenntt tt 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", distribution="elliptic",
             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, shfn=F, common=F, 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.

   distribution: If `elliptic', a multivariate elliptically-
             contoured distribution is fitted unless `pell' is
             NULL, in which case a multivariate normal distri-
             bution is fitted. If `Student  t', a multivariate
             Student t distribution is fitted and a value must
             be given for `pell'.

      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').
             If response has class, repeated, with a `linear',
             `logistic', or `pkpd' model, the covariates must
             be supplied as a Wilkinson and Rogers formula
             unless none are to be used. 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). It can-
             not be used when model is a function.

      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'. If response has class,
             repeated, with a `linear', `logistic', or `pkpd'
             model, the covariates must be supplied as a
             Wilkinson and Rogers formula unless none are to be
             used. 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.

       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
             or of the degrees of freedom parameter of the mul-
             tivariate Student t distribution. If missing and
             `distribution' is `elliptic', 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.

       shfn: If TRUE, the supplied variance function depends on
             the mean function. The name of this mean function
             must be the last argument of the variance func-
             tion.

     common: If TRUE, `mu' and `varfn' must both be functions
             with, as argument, a vector of parameters having
             some or all elements in common between them so
             that indexing is in common between them; all
             parameter estimates must be supplied in `preg'.
             If FALSE, parameters are distinct between the two
             functions and indexing starts at one in each func-
             tion.

      envir: Environment in which model formulae are to be
             interpreted or a data object of class, repeated,
             tccov, or tvcov.  If `response' has class
             `repeated', it is used as the environment.

     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.

        With two levels of nesting, the first is the individual
        and the second will consist of clusters within individ-
        uals.

        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'.)

        When an AR(1) of `exponential' form and/or a single
        random intercept is estimated for the multivariate nor-
        mal distribution, marginal and individual profiles can
        be plotted using `profile' and `iprofile' and residuals
        with `plot.residuals'.

   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', `iprofile', `kalseries', `potthoff', `profile',
        `read.list', `restovec', `rmna', `tcctomat', `tvc-
        tomat'.

   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)
        # multivariate Student t 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=5,
             distribution="Student t")
        # 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)
        # multivariate Student t 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=5,
             distribution="Student t")
        # 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)
        # multivariate Student t 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=5,
             distribution="Student t")
        #
        # 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 and Student t distributions 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)
        # finally Student t
        elliptic(y, model="logistic", preg=c(1,1,0.1,-1), pell=1,
             distribution="Student t")
        elliptic(y, model="logistic", preg=c(1,1,0.1,-1), par=0.5, pell=1,
             distribution="Student t")

