Nothing
R/gar.r
In repeated: Non-Normal Repeated Measurements Models
Defines functions
volatility.gar
volatility
print.gar
residuals.gar
fitted.gar
deviance.gar
gar
Documented in
deviance.gar
fitted.gar
gar
print.gar
residuals.gar
volatility
volatility.gar
#
# repeated : A Library of Repeated Measurements Models
# Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public Licence as published by
# the Free Software Foundation; either version 2 of the Licence, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public Licence for more details.
#
# You should have received a copy of the GNU General Public Licence
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
# SYNOPSIS
#
# gar(response=NULL, distribution="normal", times=NULL, totals=NULL,
# censor=NULL, delta=NULL, mu=NULL, shape=NULL, depend=NULL,
# shfn=FALSE, common=FALSE, preg=NULL, pshape=NULL, pdepend=NULL,
# parch=NULL, arch="square", transform="identity", link="identity",
# autocorr="exponential", order=1, envir=parent.frame(), print.level=0,
# ndigit=10, gradtol=0.00001, steptol=0.00001, fscale=1, iterlim=100,
# typsize=abs(p), stepmax=10*sqrt(p%*%p))
#
# DESCRIPTION
#
# A function to fit generalized nonlinear autoregression models with
# various distributions
##' Generalized Autoregression Models
##'
##' \code{gar} fits a first- or second-order generalized autoregression,
##' possibly with Kalman update over time (first-order only).
##'
##' Nonlinear regression models can be supplied as formulae where parameters
##' are unknowns in which case factor variables cannot be used and parameters
##' must be scalars. (See \code{\link[rmutil]{finterp}}.)
##'
##' Marginal and individual profiles can be plotted using
##' \code{\link[rmutil]{mprofile}} and \code{\link[rmutil]{iprofile}} and
##' residuals with \code{\link[rmutil]{plot.residuals}}.
##'
##' When the dispersion parameter is not constant over time, \code{volatility}
##' extracts the square root of the dispersion parameter for a fitted model.
##'
##'
##' @aliases gar volatility volatility.gar deviance.gar fitted.gar print.gar residuals.gar volatility.gar
##' @param response A list of two or three column matrices with responses,
##' corresponding times, and possibly a censor indicator, for each individual,
##' one matrix or dataframe of responses, or an object of class,
##' \code{response} (created by \code{\link[rmutil]{restovec}}) or
##' \code{repeated} (created by \code{\link[rmutil]{rmna}}). If the
##' \code{repeated} data object contains more than one response variable, give
##' that object in \code{envir} and give the name of the response variable to
##' be used here.
##' @param distribution The distribution to be fitted: binomial, Poisson,
##' exponential, negative binomial, mult Poisson, double Poisson, Consul
##' generalized Poisson, beta binomial, mult binomial, double binomial, normal,
##' inverse Gauss, logistic, gamma, Weibull, Cauchy, Laplace, Levy, Pareto,
##' beta, simplex, two-sided power, gen(eralized) gamma, gen(eralized)
##' logistic, Hjorth, Burr, gen(eralized) Weibull, gen(eralized) extreme value,
##' gen(eralized) inverse Gauss, power exponential, power variance function
##' Poisson, skew Laplace, or Student t. (For definitions of distributions, see
##' the corresponding [dpqr]distribution help.)
##' @param 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,
##' \code{response} or \code{repeated}.
##' @param totals An appropriate scalar, vector, or matrix of binomial totals
##' (only applicable for binomial, beta binomial, mult binomial, double
##' binomial). Ignored if response has class, \code{response} or
##' \code{repeated}.
##' @param censor If response is a matrix, a matrix of the same size containing
##' the censor indicator: 1=uncensored, 0=right-censored, -1=left-censored.
##' Ignored if response has class, \code{response} or \code{repeated}.
##' @param 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, \code{delta=0.01}. If the response has been
##' pretransformed, this must be multiplied by the Jacobian. This
##' transformation cannot contain unknown parameters. For example, with a log
##' transformation, \code{delta=1/y}. (The delta values for the censored
##' response are ignored.) The jacobian is calculated automatically for the
##' transform option. Ignored if response has class, \code{response} or
##' \code{repeated}.
##' @param mu A user-specified function of \code{pmu} giving the regression
##' equation for the location. It may also be 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. It must yield a value for each observation on each
##' individual.
##' @param shape An optional user-specified shape regression function; this may
##' depend on the location (function) through its second argument, in which
##' case, \code{shfn} must be set to TRUE. It may also be a formula beginning
##' with ~, specifying either a linear regression function for the shape
##' parameter in the Wilkinson and Rogers notation or a general function with
##' named unknown parameters. If it contains unknown parameters, the keyword
##' \code{mu} may be used to specify a function of the location parameter.
##' @param depend An optional user-specified regression function for the log
##' dependence parameter. It may also be a formula beginning with ~, specifying
##' either a linear regression function for the dependence parameter in the
##' Wilkinson and Rogers notation or a general function with named unknown
##' parameters. If used, \code{order} must be one.
##' @param shfn If TRUE, the supplied shape function depends on the location
##' function. The name of this location function must be the last argument of
##' the shape function.
##' @param common If TRUE, \code{mu} and \code{shape} must both be either
##' 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
##' or formulae with unknowns. All parameter estimates must be supplied in
##' \code{preg}. If FALSE, parameters are distinct between the two functions
##' and indexing starts at one in each function.
##' @param preg The initial parameter estimates for the location regression
##' function. If \code{mu} is a formula with unknown parameters, their
##' estimates must be supplied either in their order of appearance in the
##' expression or in a named list.
##' @param pdepend One or two estimates of the dependence parameters for the
##' Kalman update. With one, it is Markovian and, with two, it is
##' nonstationary. For the latter, the \code{order} must be one. If
##' \code{depend} is a function or formula, the corresponding number of
##' estimates must be supplied. Either pdepend or parch or both must be
##' supplied.
##' @param parch Estimate for an ARCH model where the shape parameter depends
##' on the square of the previous residual. Either pdepend or parch or both
##' must be supplied.
##' @param arch If \code{square}, then \code{shape+parch^diff*residual^2}; if
##' \code{absolute value}, then \code{shape+parch^diff*|residual|}; if
##' \code{exponential}, then \code{shape*exp(parch*residual^2*diff)}, where
##' \code{diff} is the length of time since the previous observation and
##' \code{residual} is the previous residual or innovation.
##' @param pshape Zero to two estimates for the shape parameters, depending on
##' the distribution, if \code{shape} is not a function; otherwise, estimates
##' for the parameters in this function, with one extra at the end for
##' three-parameter distributions. If \code{shape} is a formula with unknown
##' parameters, their estimates must be supplied either in their order of
##' appearance in the expression or in a named list (only for two-parameter
##' distributions).
##' @param transform Transformation of the response variable: \code{identity},
##' \code{exp}, \code{square}, \code{sqrt}, or \code{log}.
##' @param link Link function for the mean: \code{identity}, \code{exp},
##' \code{square}, \code{sqrt}, \code{log}, \code{logit}, \code{cloglog} or
##' \code{loglog} (last three only for binary/binomial-type data).
##' @param autocorr The form of the (second if two) dependence function:
##' \code{exponential} is the usual \eqn{\rho^{|t_i-t_j|}}{rho^|t_i-t_j|};
##' \code{gaussian} is \eqn{\rho^{(t_i-t_j)^2}}{rho^((t_i-t_j)^2)};
##' \code{cauchy} is \eqn{1/(1+\rho(t_i-t_j)^2)}{1/(1+rho(t_i-t_j)^2)};
##' \code{spherical} is
##' \eqn{((|t_i-t_j|\rho)^3-3|t_i-t_j|\rho+2)/2}{((|t_i-t_j|rho)^3-3|t_i-t_j|rho+2)/2}
##' for \eqn{|t_i-t_j|\leq1/\rho}{|t_i-t_j|<=1/rho} and zero otherwise;
##' \code{IOU} is the integrated Ornstein-Uhlenbeck process, \eqn{(2\rho
##' \min(t_i,t_j)+\exp(-\rho t_i) }{(2rho min(t_i,t_j)+exp(-rho t_i) +exp(-rho
##' t_j)-1 -exp(rho|ti-t_j|))/2rho^3}\eqn{+\exp(-\rho t_j)-1
##' -\exp(\rho|ti-t_j|))/2\rho^3}{(2rho min(t_i,t_j)+exp(-rho t_i) +exp(-rho
##' t_j)-1 -exp(rho|ti-t_j|))/2rho^3}.
##' @param order First- or second-order stationary autoregression.
##' @param envir Environment in which model formulae are to be interpreted or a
##' data object of class, \code{repeated}, \code{tccov}, or \code{tvcov}; the
##' name of the response variable should be given in \code{response}. If
##' \code{response} has class \code{repeated}, it is used as the environment.
##' @param print.level Arguments for nlm.
##' @param typsize Arguments for nlm.
##' @param ndigit Arguments for nlm.
##' @param gradtol Arguments for nlm.
##' @param stepmax Arguments for nlm.
##' @param steptol Arguments for nlm.
##' @param iterlim Arguments for nlm.
##' @param fscale Arguments for nlm.
##' @return A list of classes \code{gar} and \code{recursive} is returned that
##' contains all of the relevant information calculated, including error codes.
##'
##' The volatility vector for models with a shape regression function and ARCH
##' models contains the square root of the dispersion parameter at each time
##' point.
##' @author J.K. Lindsey
### @seealso \code{\link[growth]{carma}}, \code{\link[growth]{elliptic}},
### \code{\link[rmutil]{finterp}}, \code{\link[repeated]{gnlmm}},
### \code{\link[gnlm]{gnlr}}, \code{\link[rmutil]{iprofile}},
### \code{\link[repeated]{kalcount}}, \code{\link[repeated]{kalseries}},
### \code{\link[event]{kalsurv}}, \code{\link[rmutil]{mprofile}},
### \code{\link[rmutil]{plot.residuals}}, \code{\link[rmutil]{read.list}},
### \code{\link[rmutil]{restovec}}, \code{\link[rmutil]{rmna}},
### \code{\link[rmutil]{tcctomat}}, \code{\link[rmutil]{tvctomat}}.
##' @references Lindsey, J.K. (1997) Applying Generalized Linear Models.
##' Springer, pp.\ 93--101
##'
##' Lambert, P. (1996) Statistics in Medicine 15, 1695-1708
##' @keywords models
##' @examples
##'
##' # first-order one-compartment model
##' # data objects for formulae
##' dose <- c(2,5)
##' dd <- tcctomat(dose)
##' times <- matrix(rep(1:20,2), nrow=2, byrow=TRUE)
##' tt <- tvctomat(times)
##' # vector covariates for functions
##' dose <- c(rep(2,20),rep(5,20))
##' times <- rep(1:20,2)
##' # functions
##' 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)
##' # response
##' conc <- matrix(rgamma(40,shape(log(c(0.1,0.4))),
##' scale=mu(log(c(1,0.3,0.2))))/shape(log(c(0.1,0.4))),ncol=20,byrow=TRUE)
##' conc[,2:20] <- conc[,2:20]+0.5*(conc[,1:19]-matrix(mu(log(c(1,0.3,0.2))),
##' ncol=20,byrow=TRUE)[,1:19])
##' conc <- restovec(ifelse(conc>0,conc,0.01),name="conc")
##' reps <- rmna(conc, ccov=dd, tvcov=tt)
##' # constant shape parameter
##' gar(conc, dist="gamma", times=1:20, mu=mu,
##' preg=log(c(1,0.4,0.1)), pdepend=0.5, pshape=1)
##' \dontrun{# or
##' gar(conc, dist="gamma", times=1:20, mu=~exp(absorption-volume)*
##' dose/(exp(absorption)-exp(elimination))*
##' (exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
##' preg=list(absorption=1,elimination=log(0.4),volume=log(0.1)),
##' pdepend=0.5, pshape=1, envir=reps)
##' # generalized gamma distribution
##' gar(conc, dist="gen gamma", times=1:20, mu=mu,
##' preg=log(c(1,0.4,0.1)), pdepend=0.3, pshape=c(.1,1))
##' # (if the covariates contained NAs, reps would have to be used as
##' # response instead of conc)
##' #
##' # time dependent shape parameter
##' gar(conc, dist="gamma", times=1:20, mu=mu, shape=shape,
##' preg=log(c(1,0.4,0.1)), pdepend=0.25, pshape=c(exp(-2),exp(-.57)))
##' # or
##' gar(conc, dist="gamma", times=1:20, mu=~exp(absorption-volume)*
##' dose/(exp(absorption)-exp(elimination))*
##' (exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
##' shape=~exp(b1-b2)*times*dose*exp(-exp(b1)*times),
##' preg=list(absorption=0,elimination=log(0.4),volume=log(0.1)),
##' pdepend=0.3, pshape=list(b1=exp(-2),b2=exp(-.57)), envir=reps)
##' # generalized gamma distribution
##' gar(conc, dist="gen gamma", times=1:20, mu=mu, shape=shape,
##' preg=log(c(1,0.4,0.1)), pdepend=0.5,
##' pshape=c(exp(-2),exp(-.57),2))
##' #
##' # shape function depends on location parameter
##' shape <- function(p, mu) p[1]+p[2]*mu
##' gar(conc, dist="gamma", times=1:20, mu=mu, shape=shape, shfn=TRUE,
##' preg=log(c(1,0.4,.10)), pdepend=0.15, pshape=c(1,2))
##' # or
##' gar(conc, dist="gamma", times=1:20, mu=mu, shape=~a+d*mu, shfn=TRUE,
##' preg=log(c(1,0.4,.10)), pdepend=0.15, pshape=c(1,2))
##'
##' }
##' @export gar
gar <- function(response=NULL, distribution="normal", times=NULL, totals=NULL,
censor=NULL, delta=NULL, mu=NULL, shape=NULL, depend=NULL, shfn=FALSE,
common=FALSE, preg=NULL, pshape=NULL, pdepend=NULL, parch=NULL,
arch="square", transform="identity", link="identity",
autocorr="exponential", order=1, envir=parent.frame(), print.level=0,
ndigit=10, gradtol=0.00001, steptol=0.00001, fscale=1, iterlim=100,
typsize=abs(p), stepmax=10*sqrt(p%*%p)){
#
# set up likelihood to call C code
#
likekal <- function(p){
eta <- mu1(p)
if(sh){
shr <- sh1(p)
theta <- exp(p[np])
if(!tm)theta <- c(p[npr1:nprd],theta)}
else theta <- if(tm)p[nprd1:np] else p[npr1:np]
if(tm)dep <- dep1(p)
if(archtype>0)parch <- p[npa]
z <- .C("gar_c",
y=as.double(y),
total=as.double(response$response$n),
my=as.integer(3*max(y)),
nobs=as.integer(nobs(response)),
nind=as.integer(nind),
times=as.double(response$response$times),
censor=as.integer(censor),
cens=as.integer(!is.null(response$response$censor)),
eta=as.double(eta),
tm=as.integer(tm),
theta=as.double(theta),
dep=as.double(dep),
model=as.integer(mdl),
thp=as.integer(thp),
shape=as.double(shr),
sh=as.integer(sh),
link=as.integer(lnk),
ar=as.integer(ar),
order=as.integer(order),
parch=as.double(parch),
arch=as.integer(archtype),
pred=double(n),
rpred=double(n),
volatility=double(n),
like=double(1),
#DUP=FALSE,
PACKAGE="repeated")
z$like+jacob}
call <- sys.call()
#
# check model required
#
tmp <- c("binomial","Poisson","exponential","negative binomial",
"mult Poisson","double Poisson","Consul","beta binomial",
"mult binomial","double binomial","normal","logistic","Cauchy",
"Weibull","gamma","Laplace","inverse Gauss","Pareto","Levy",
"beta","simplex","two-sided power","gen gamma","gen logistic","Hjorth",
"Burr","gen Weibull","gen extreme value","gen inverse Gauss",
"power exponential","power variance function Poisson",
"skew Laplace","Student t")
mdl <- match(distribution <- match.arg(distribution,tmp),tmp)
bmdl <- mdl==1||mdl==8||mdl==9||mdl==10
par1 <- distribution=="binomial"||distribution=="Poisson"||
distribution=="exponential"
par3 <- mdl>22
tmp <- c("exponential","gaussian","cauchy","spherical","IOU")
ar <- match(autocorr <- match.arg(autocorr,tmp),tmp)
#
# check transformation and link function
#
transform <- match.arg(transform,c("identity","exp","square","sqrt","log"))
if(transform!="identity"&&((mdl<11&&mdl!=3)||mdl==31))
stop("transformations not valid for discrete data models")
tmp <- c("identity","exp","square","sqrt","log","logit","cloglog","loglog")
lnk <- match(link <- match.arg(link,tmp),tmp)
if((link=="logit"||link=="cloglog"||link=="loglog")&&!bmdl)
stop("logit, cloglog, and loglog links can only be used with binary data")
#
# check that either AR or ARCH is being fitted
#
if(is.null(pdepend)&&is.null(parch))
stop("Either pdepend or parch must be supplied")
if(!is.null(parch)){
if(par1)stop(paste("arch model not possible for",distribution,"distribution"))
if(length(parch)!=1)stop("parch should contain one value")
tmp <- c("square","absolute value","exponential")
archtype <- match(match.arg(arch,tmp),tmp)
if(archtype!=3&&parch<=0)
stop("parch must be positive for square and absolute value models")}
else archtype <- 0
#
# check if a data object is being supplied
#
type <- "unknown"
respenv <- exists(deparse(substitute(response)),envir=parent.frame())&&
inherits(response,"repeated")
if(respenv){
if(dim(response$response$y)[2]>1)
stop("gar only handles univariate responses")
if(!is.null(response$NAs)&&any(response$NAs))
stop("gar does not handle data with NAs")
type <- response$response$type}
envname <- if(respenv)deparse(substitute(response))
else if(!is.null(class(envir)))deparse(substitute(envir))
else NULL
#
#if response is not a data object, find or make one
#
if(!respenv){
if(inherits(envir,"repeated")){
# if envir, remove extra (multivariate) responses
if(!is.null(envir$NAs)&&any(envir$NAs))
stop("gar does not handle data with NAs")
cn <- deparse(substitute(response))
if(length(grep("\"",cn))>0)cn <- response
if(length(cn)>1)stop("only one response variable allowed")
response <- envir
col <- match(cn,colnames(response$response$y))
if(is.na(col))stop(paste("response variable",cn,"not found"))
type <- response$response$type[col]
if(dim(response$response$y)[2]>1){
response$response$y <- response$response$y[,col,drop=FALSE]
if(!is.null(response$response$delta)){
response$response$delta <- response$response$delta[,col,drop=FALSE]
if(all(response$response$delta==1)||all(is.na(response$response$delta)))response$response$delta <- NULL}
if(!is.null(response$response$n)){
response$response$n <- response$response$n[,col,drop=FALSE]
if(all(is.na(response$response$n)))response$response$n <- NULL}
if(!is.null(envir$response$censor)){
response$response$censor <- response$response$censor[,col,drop=FALSE]
if(all(is.na(response$response$censor)))response$response$censor <- NULL}}}
else {
if(!inherits(response,"response")){
if(is.null(times)){
if(is.matrix(response))times <- matrix(1:dim(response)[2],nrow=dim(response)[1],ncol=dim(response)[2],byrow=T)
else if(is.vector(response))times <- 1:length(response)}
response <- if(bmdl)
restovec(response, times=times, totals=totals)
else restovec(response, times=times, censor=censor, delta=delta)}
type <- response$type
response <- rmna(response=response)}}
if(is.null(times(response)))stop("No times available")
if((inherits(envir,"repeated")&&
(length(nobs(response))!=length(nobs(envir))||
any(nobs(response)!=nobs(envir))))||(inherits(envir,"tvcov")&&
(length(nobs(response))!=length(envir$tvcov$nobs)||
any(nobs(response)!=envir$tvcov$nobs))))
stop("response and envir objects are incompatible")
if(bmdl){
if(type!="unknown"&&type!="nominal")stop("nominal data required")
if(is.null(response$response$n)){
if(any(response$response$y!=0&&response$response$y!=1))
stop("responses must be binary if totals are not supplied")
else response$response$n <- rep(1,dim(response$response$y)[1])}}
y <- response$response$y
n <- dim(y)[1]
#
# if a data object was supplied, modify formulae or functions to read from it
#
mu3 <- sh3 <- dep3 <- dep <- NULL
if(respenv||inherits(envir,"repeated")||inherits(envir,"tccov")||inherits(envir,"tvcov")){
if(is.null(envname))envname <- deparse(substitute(envir))
# location function or formula
if(inherits(mu,"formula")){
mu3 <- if(respenv)finterp(mu,.envir=response,.name=envname)
else finterp(mu,.envir=envir,.name=envname)}
else if(is.function(mu)){
if(is.null(attr(mu,"model"))){
tmp <- parse(text=deparse(mu)[-1])
mu <- if(respenv)fnenvir(mu,.envir=response,.name=envname)
else fnenvir(mu,.envir=envir,.name=envname)
mu3 <- mu
attr(mu3,"model") <- tmp}
else mu3 <- mu}
# shape function or formula
if(inherits(shape,"formula")){
sh3 <- if(respenv)finterp(shape,.envir=response,.name=envname)
else finterp(shape,.envir=envir,.name=envname)}
else if(is.function(shape)){
if(is.null(attr(shape,"model"))){
tmp <- parse(text=deparse(shape)[-1])
shape <- if(respenv)fnenvir(shape,.envir=response,.name=envname)
else fnenvir(shape,.envir=envir,.name=envname)
sh3 <- shape
if(shfn){
pos <- grep("mu",attr(sh3,"parameters"))
attr(sh3,"parameters") <- attr(sh3,"parameters")[-pos]}
attr(sh3,"model") <- tmp}
else sh3 <- shape}
# AR function or formula
if(inherits(depend,"formula")){
dep3 <- if(respenv)finterp(depend,.envir=response,.name=envname)
else finterp(depend,.envir=envir,.name=envname)}
else if(is.function(depend)){
if(is.null(attr(depend,"model"))){
tmp <- parse(text=deparse(depend)[-1])
depend <- if(respenv)fnenvir(depend,.envir=response,.name=envname)
else fnenvir(depend,.envir=envir,.name=envname)
dep3 <- depend
attr(dep3,"model") <- tmp}
else dep3 <- depend}}
#
# transform location formula to function and check number of parameters
#
npr <- length(preg)
npr1a <- npr1 <- npr+1
if(inherits(mu,"formula")){
mu2 <- if(respenv)finterp(mu,.envir=response,.name=envname)
else finterp(mu,.envir=envir,.name=envname)
npt1 <- length(attr(mu2,"parameters"))
if(is.character(attr(mu2,"model"))){
# W&R formula
if(length(attr(mu2,"model"))==1){
# intercept model
mu1 <- function(p) p[1]*rep(1,n)
attributes(mu1) <- attributes(mu2)
mu2 <- NULL}}
else {
# formula with unknowns
if(npr!=npt1&&!common){
cat("\nParameters are ")
cat(attr(mu2,"parameters"),"\n")
stop(paste("preg should have",npt1,"estimates"))}
if(is.list(preg)){
if(!is.null(names(preg))){
o <- match(attr(mu2,"parameters"),names(preg))
preg <- unlist(preg)[o]
if(sum(!is.na(o))!=length(preg))stop("invalid estimates for mu - probably wrong names")}
else preg <- unlist(preg)}}
if(!is.null(mu2)){
if(inherits(envir,"tccov")){
cv <- covind(response)
mu1 <- function(p) mu2(p)[cv]
attributes(mu1) <- attributes(mu2)}
else {
mu1 <- mu2
rm(mu2)}}}
else if(is.null(mu)){
mu1 <- function(p) p[1]*rep(1,n)
npt1 <- 1}
else mu1 <- mu
#
# give appropriate attributes to mu1 for printing
#
if(!is.null(mu1)&&is.null(attr(mu1,"parameters"))){
attributes(mu1) <- if(is.function(mu)){
if(!inherits(mu,"formulafn")){
if(respenv)attributes(fnenvir(mu,.envir=response))
else attributes(fnenvir(mu,.envir=envir))}
else attributes(mu)}
else {
if(respenv)attributes(fnenvir(mu1,.envir=response))
else attributes(fnenvir(mu1,.envir=envir))}}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(mu1))length(attr(mu1,"parameters"))
else if(is.null(mu1))NULL
else npt1
if(!is.null(nlp)&&!common&&nlp!=npr)
stop(paste("preg should have",nlp,"initial estimates"))
#
# transform shape formula to function and check number of parameters
#
nprd <- npr+length(pdepend)
nprd1 <- if(common&&!inherits(shape,"formula")) 1 else nprd+1
nps1 <- nps <- length(pshape)
if(par3)nps1 <- nps1-1
if(inherits(shape,"formula")){
old <- if(common)mu1 else NULL
mufn <- if(shfn)"mu" else NULL
sh3 <- if(respenv)finterp(shape,.envir=response,.start=nprd1,.name=envname,.old=old,.args=mufn)
else finterp(shape,.envir=envir,.start=nprd1,.name=envname,.old=old,.args=mufn)
tmp <- attributes(sh3)
sh2 <- if(shfn)function(p) sh3(p,mu1(p)) else sh3
attributes(sh2) <- tmp
npt2 <- length(attr(sh2,"parameters"))
if(is.character(attr(sh2,"model"))){
# W&R formula
if(length(attr(sh2,"model"))==1){
# intercept model
sh1 <- function(p) p[nprd1]*rep(1,n)
attributes(sh1) <- attributes(sh2)
sh2 <- NULL}}
else {
# formula with unknowns
if(nps1!=npt2&&!common){
cat("\nParameters are ")
cat(attr(sh2,"parameters"),"\n")
stop(paste("pshape should have",npt2+par3,"estimates"))}
if(is.list(pshape)){
if(!is.null(names(pshape))){
o <- match(attr(sh2,"parameters"),names(pshape))
pshape <- unlist(pshape)[o]
if(sum(!is.na(o))!=length(pshape))stop("invalid estimates for shape - probably wrong names")}
else pshape <- unlist(pshape)}}
if(!is.null(sh2)){
if(inherits(envir,"tccov")){
cv <- covind(response)
sh1 <- function(p) sh2(p)[cv]
attributes(sh1) <- attributes(sh2)}
else {
sh1 <- sh2
rm(sh2)}}}
else if(!is.function(shape)&&!par1){
sh1 <- function(p) p[nprd1]*rep(1,n)
npt2 <- 1
if(length(pshape)!=1+par3)
stop(paste("pshape must contain",1+par3,"estimate(s)"))}
else sh1 <- if(shfn)function(p) shape(p[nprd1:np], mu1(p))
else function(p) shape(p[nprd1:np])
#
# give appropriate attributes to sh1 for printing
#
if(!is.null(sh1)&&is.null(attr(sh1,"parameters")))
attributes(sh1) <- if(is.function(shape)){
if(!inherits(shape,"formulafn")){
if(respenv)attributes(fnenvir(shape,.envir=response))
else attributes(fnenvir(shape,.envir=envir))}
else attributes(shape)}
else {
if(respenv)attributes(fnenvir(sh1,.envir=response))
else attributes(fnenvir(sh1,.envir=envir))}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(shape))length(attr(sh1,"parameters"))-shfn
else if(is.null(shape))NULL
else npt2
if(!is.null(nlp)&&!common&&nlp!=nps1)
stop(paste("pshape should have",nlp+par3,"initial estimates"))
if(common){
nlp <- length(unique(c(attr(mu1,"parameters"),attr(sh1,"parameters"))))
if(nlp!=npr)stop(paste("with a common parameter model, preg should contain",nlp,"estimates"))}
#
# transform AR formula to function and check number of parameters
#
npd <- length(pdepend)
if(inherits(depend,"formula")){
dep2 <- if(respenv)finterp(depend,.envir=response,.start=npr1,.name=envname)
else finterp(depend,.envir=envir,.start=npr1,.name=envname)
npt3 <- length(attr(dep2,"parameters"))
if(is.character(attr(dep2,"model"))){
# W&R formula
if(length(attr(dep2,"model"))==1){
# intercept model
dep1 <- function(p) p[npr1]*rep(1,n)
attributes(dep1) <- attributes(dep2)
dep2 <- NULL}}
else {
# formula with unknowns
if(npd!=npt3){
cat("\nParameters are ")
cat(attr(dep2,"parameters"),"\n")
stop(paste("pdepend should have",npt3,"estimates"))}
if(is.list(pdepend)){
if(!is.null(names(pdepend))){
o <- match(attr(dep2,"parameters"),names(pdepend))
pdepend <- unlist(pdepend)[o]
if(sum(!is.na(o))!=length(pdepend))stop("invalid estimates for depend - probably wrong names")}
else pdepend <- unlist(pdepend)}}
if(!is.null(dep2)){
if(inherits(envir,"tccov")){
cv <- covind(response)
dep1 <- function(p) dep2(p)[cv]
attributes(dep1) <- attributes(dep2)}
else {
dep1 <- dep2
rm(dep2)}}}
else if(is.function(dep))dep1 <- dep
else {
dep1 <- NULL
npt3 <- 1}
#
# give appropriate attributes to dep1 for printing
#
if(!is.null(dep1)&&is.null(attr(dep1,"parameters")))
attributes(dep1) <- if(is.function(depend)){
if(!inherits(depend,"formulafn")){
if(respenv)attributes(fnenvir(depend,.envir=response))
else attributes(fnenvir(depend,.envir=envir))}
else attributes(depend)}
else {
if(respenv)attributes(fnenvir(dep1,.envir=response))
else attributes(fnenvir(dep1,.envir=envir))}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(depend))length(attr(dep1,"parameters"))
else if(is.null(depend))NULL
else npt3
if(!is.null(nlp)&&nlp!=npd)
stop(paste("pdepend should have",nlp,"initial estimates"))
#
# check response and functions
#
if(mdl<11){
if((distribution=="Poisson"||distribution=="negative binomial"
||distribution=="mult Poisson"
||distribution=="double Poisson"||distribution=="Consul")
&&type!="unknown"&&type!="discrete")
stop("discrete data required")
if(distribution=="exponential"&&type!="unknown"&&type!="duration")
stop("duration data required")
if(any(y<0))stop("All responses must be non-negative")}
else if(mdl==20||mdl==21||mdl==22){
if(type!="unknown"&&type!="continuous")stop("continuous data required")
if(any(y<=0)||any(y>=1))stop("All responses must be between 0 and 1")}
else if(mdl==31){
if(type!="unknown"&&type!="discrete")stop("discrete data required")
if(any(y<0))stop("All responses must be non-negative")}
else if(mdl!=11&&mdl!=12&&mdl!=13&&mdl!=16&&mdl!=19&&mdl!=24&&mdl!=30&&
mdl!=32&&mdl!=33){
if(type!="unknown"&&type!="duration"&&type!="continuous")
stop("duration data required")
if(any(y<=0))stop("All responses must be positive")}
else {
if(type!="unknown"&&type!="continuous"&&type!="duration")
stop("continuous data required")
if(distribution=="Levy"&&any(response$response$y<=mu1(preg)))
stop("Location function must give values strictly less than corresponding observations")}
if(distribution=="Pareto"&&any(sh1(pshape)<=1))
stop("shape function must be > 1")
#
# set up censoring indicator
#
censor <- response$response$censor
if(is.null(censor))censor <- rep(1,n)
else if(bmdl||mdl==2||mdl==4||mdl==5||mdl==6||mdl==7||mdl==31)
stop(paste("Censored data not allowed for",distribution,"distribution"))
#
# set up transformation and Jacobian
#
nind <- length(nobs(response))
if(transform=="identity")jacob <- 0
else if(transform=="exp"){
jacob <- -sum(y[censor==1])
y <- exp(y)}
else if(any(y==0))stop("Zero response values: invalid transformation")
else if(transform=="square"){
jacob <- -sum(log(abs(y[censor==1])))-n*log(2)
y <- y^2}
else if(any(y<0))stop("Nonpositive response values: invalid transformation")
else if(transform=="sqrt"){
jacob <- sum(log(y[censor==1]))/2+n*log(2)
y <- sqrt(y)}
else if(transform=="log"){
jacob <- sum(log(y[censor==1]))
y <- log(y)}
#
# set up unit of measurement
#
if(!is.null(response$response$delta)){
if(length(response$response$delta)==1)
jacob <- jacob-length(y[censor==1])*log(response$response$delta)
else jacob <- jacob-sum(log(response$response$delta[censor==1]))}
#
# check AR model
#
if(is.null(pdepend))order <- 0
else if(order!=1&&order!=2)stop("Autoregression must have order 1 or 2")
tm <- !is.null(depend)
if(tm)order <- 0
if(order==2&&length(pdepend)!=2)
stop("2 estimates of dependence parameters must be given")
else if(order>0&&!tm&&length(pdepend)!=1&&length(pdepend)!=2)
stop("One or two estimates of dependence parameters must be given")
thp <- length(pdepend)==2&&order==1&&!tm
if(!tm&&any(pdepend<=0))stop("All dependence parameters must be positive")
#
# check location function
#
if(length(mu1(preg))!=dim(response$response$y)[1])
stop("The mu function must provide an estimate for each observation")
else if(any(is.na(mu1(preg))))
stop("Non-numerical mu: probably invalid initial values")
#
# set up parameter vector
#
p <- if(tm)c(preg,pdepend) else if(order>0)c(preg,-log(pdepend)) else preg
sh <- mdl>3&&(is.function(shape)||inherits(shape,"formula"))
if(!par1){
if(!sh)p <- c(p,if(mdl==31)c(log(pshape[1]),pshape[2])else log(pshape))
else {
if(par3){
if(pshape[nps]<=0&&mdl!=31)
stop("extra shape parameter must be positive")
p <- if(common)c(p,if(mdl==31)pshape[nps]
else log(pshape[nps]))
else c(p,pshape[1:(nps-1)],if(mdl==31)
pshape[nps] else log(pshape[nps]))}
else p <- c(p,pshape)}}
np <- length(p)
if(archtype>0)p <- c(p,if(archtype==1)log(parch)else parch)
npa <- length(p)
#
# check shape function
#
if(!sh){
if((mdl<=3&&nps!=0)||(mdl>3&&mdl<23&&nps!=1)||(mdl>=23&&nps!=2))
stop("Incorrect number of shape parameter estimates")
else if(nps>0){
if(mdl!=31&&any(pshape<=0)||mdl==31&&pshape[1]<=0)
stop("All shape parameters must be positive")}
shr <- rep(0,dim(response$response$y)[1])}
else {
if(any(is.na(sh1(p))))
stop("The shape model returns NAs: probably invalid initial values")
if(length(sh1(p))!=dim(response$response$y)[1])
stop("The shape function must provide an estimate for each observation")}
#
# check likelihood function and optimize with nlm
#
if(fscale==1)fscale <- likekal(p)
if(is.na(likekal(p)))
stop("Likelihood returns NAs: probably invalid initial values")
z0 <- nlm(likekal, p, hessian=TRUE, print.level=print.level,
typsize=typsize, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
steptol=steptol, iterlim=iterlim, fscale=fscale)
like <- z0$minimum
#
# calculate se's
#
a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
else qr(z0$hessian)$rank
if(a==npa)cov <- solve(z0$hessian)
else cov <- matrix(NA,ncol=npa,nrow=npa)
se <- sqrt(diag(cov))
corr <- cov/(se%o%se)
dimnames(corr) <- list(1:npa,1:npa)
#
# back transform parameters
#
eta <- mu1(z0$estimate)
if(sh){
shr <- sh1(z0$estimate)
theta <- exp(z0$estimate[np])
if(!tm)theta <- c(z0$estimate[npr1:nprd],theta)}
else theta <- if(tm)z0$estimate[nprd1:np] else z0$estimate[npr1:np]
if(tm)dep <- dep1(z0$estimate)
#
# calculate predicted values (and volatility)
#
if(archtype>0)parch <- z0$estmate[npa]
z <- .C("gar_c",
y=as.double(y),
total=as.double(response$response$n),
my=as.integer(3*max(y)),
nobs=as.integer(nobs(response)),
nind=as.integer(nind),
times=as.double(response$response$times),
censor=as.integer(censor),
cens=as.integer(!is.null(response$response$censor)),
eta=as.double(eta),
tm=as.integer(tm),
theta=as.double(theta),
dep=as.double(dep),
model=as.integer(mdl),
thp=as.integer(thp),
shape=as.double(shr),
sh=as.integer(sh),
link=as.integer(lnk),
ar=as.integer(ar),
order=as.integer(order),
parch=as.double(parch),
arch=as.integer(archtype),
pred=double(n),
rpred=double(n),
volatility=double(n),
like=double(1),
#DUP=FALSE,
PACKAGE="repeated")
if(archtype==0)z$volatility <- NULL
#
# return object
#
if(!is.null(mu3))mu1 <- mu3
if(!is.null(sh3))sh1 <- sh3
if(!is.null(dep3))dep1 <- dep3
z <- list(
call=call,
distribution=distribution,
mu=mu1,
formula=mu,
shape=shape,
depend=depend,
sh1=sh1,
shfn=shfn,
common=common,
dep1=dep1,
response=response$response,
link=link,
order=order,
autocorr=autocorr,
transform=transform,
maxlike=like,
aic=like+npa,
df=dim(response$response$y)[1]-npa,
np=np,
npa=npa,
npr=npr,
nps=nps,
npd=npd,
thp=thp,
tm=tm,
archtype=archtype,
coefficients=z0$estimate,
se=se,
cov=cov,
corr=corr,
pred=z$pred,
rpred=z$rpred,
volatility=z$volatility,
grad=z0$gradient,
iterations=z0$iterations,
code=z0$code)
class(z) <- c("gar","recursive")
return(z)}
### standard methods
###
#' @export
deviance.gar <- function(object, ...) 2*object$maxlike
#' @export
fitted.gar <- function(object, recursive=TRUE, ...) if(recursive) object$rpred else object$pred
#' @export
residuals.gar <- function(object, recursive=TRUE, ...){
z <- object
if(z$transform=="exp")z$response$y <- exp(z$response$y)
else if(z$transform=="square")z$response$y <- z$response$y^2
else if(z$transform=="sqrt")z$response$y <- sqrt(z$response$y)
else if(z$transform=="log")z$response$y <- log(z$response$y)
if(recursive) z$response$y-z$rpred else z$response$y-z$pred}
### print method
###
#' @export
print.gar <- function(x,digits=max(3,.Options$digits-3),correlation=TRUE, ...){
z <- x
np1 <- if(z$distribution=="binomial"||z$distribution=="exponential"
||z$distribution=="Poisson") 0
else if(z$distribution=="gen gamma"
||z$distribution=="gen logistic"
||z$distribution=="Hjorth"||z$distribution=="Burr"
||z$distribution=="gen Weibull"
||z$distribution=="gen extreme value"
||z$distribution=="gen inverse Gauss"
||z$distribution=="power exponential"
||z$distribution=="power variance function Poisson"
||z$distribution=="skew Laplace"
||z$distribution=="Student t") 2
else 1
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
cat("Number of subjects ",length(nobs(z)),"\n")
cat("Number of observations",length(z$response$y),"\n")
cat("Transformation ",z$trans,"\n")
cat("Link function ",z$link,"\n\n")
cat(z$distribution,"distribution\n")
if(z$order<2)cat("First order ")
else cat("Second order ")
cat(z$autocorr,"dependence\n")
cat("\n-Log likelihood ",z$maxlike,"\n")
cat("Degrees of freedom",z$df,"\n")
cat("AIC ",z$aic,"\n")
cat("Iterations ",z$iterations,"\n\n")
if(z$common)cat("Location model\n")
else cat("Location parameters\n")
if(!is.null(attr(z$mu,"formula")))
cat(deparse(attr(z$mu,"formula")),sep="\n")
else if(!is.null(attr(z$mu,"model"))){
t <- deparse(attr(z$mu,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- if(is.character(attr(z$mu,"model")))attr(z$mu,"model")
else attr(z$mu,"parameters")
coef.table <- cbind(z$coef[1:z$npr],z$se[1:z$npr])
colname <- c("estimate","se")
if(!z$common){
dimnames(coef.table) <- list(cname,colname)
print.default(coef.table, digits=digits, print.gap=2)}
else {
cat("\nShape model\n")
if(!is.null(attr(z$sh1,"formula")))
cat(deparse(attr(z$sh1,"formula")),sep="\n")
else if(!is.null(attr(z$sh1,"model"))){
t <- deparse(attr(z$sh1,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- c(cname,if(is.character(attr(z$sh1,"model")))
attr(z$sh1,"model")
else attr(z$sh1,"parameters")[1:(length(attr(z$sh1,"parameters"))-z$shfn)])
cname <- unique(cname)
cat("\nCommon parameters\n")
dimnames(coef.table) <- list(cname,colname)
print.default(coef.table, digits=digits, print.gap=2)}
if(z$order>0||z$tm){
if(z$thp||z$order==2){
cat("\nDependence parameters\n")
if(z$thp)cname <- c("phi","rho")
else cname <- c("rho1","rho2")
coef.table <- cbind(z$coef[(z$npr+1):(z$npr+2)],
z$se[(z$npr+1):(z$npr+2)],
exp(-z$coef[(z$npr+1):(z$npr+2)]))}
else if(z$tm){
cat("\nDependence parameters\n")
if(!is.null(attr(z$dep1,"formula")))
cat(deparse(attr(z$dep1,"formula")),sep="\n")
else if(!is.null(attr(z$dep1,"model"))){
t <- deparse(attr(z$dep1,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- if(is.character(attr(z$dep1,"model")))
attr(z$dep1,"model")
else attr(z$dep1,"parameters")
coef.table <- cbind(z$coef[(z$npr+1):(z$npr+z$npd)],
z$se[(z$npr+1):(z$npr+z$npd)])
colname <- c("estimate","se")
dimnames(coef.table) <- list(cname,colname)}
else {
cat("\nDependence parameter\n")
cname <- "rho"
coef.table <- cbind(z$coef[z$npr+1],
z$se[z$npr+1],exp(-z$coef[z$npr+1]))}
if(!z$tm)dimnames(coef.table) <- list(cname, c("estimate","se","parameter"))
print.default(coef.table, digits=digits, print.gap=2)}
if(np1>0&&!z$common){
cat("\nShape parameters\n")
if(is.null(z$shape)){
cname <- "shape"
if(np1==2)cname <- c(cname,"psi")
coef.table <- cbind(z$coef[(z$np-np1+1):z$np],
z$se[(z$np-np1+1):z$np],if(z$distribution==
"power variance function Poisson")
c(exp(z$coef[(z$np-np1+1):(z$np-1)]),z$coef[z$np]) else
exp(z$coef[(z$np-np1+1):z$np]))
dimnames(coef.table) <- list(cname, c("estimate","se","parameter"))}
else {
if(!is.null(attr(z$sh1,"formula")))
cat(deparse(attr(z$sh1,"formula")),sep="\n")
else if(!is.null(attr(z$sh1,"model"))){
t <- deparse(attr(z$sh1,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- if(is.character(attr(z$sh1,"model")))
attr(z$sh1,"model")
else attr(z$sh1,"parameters")
np2 <- length(cname)+(np1==2)
coef.table <- cbind(z$coef[(z$np-np2+1):z$np],
z$se[(z$np-np2+1):z$np])
if(np1==2){
cname <- c(cname,"psi")
coef.table <- cbind(coef.table,
c(coef.table[1:(dim(coef.table)[1]-1),1],
if(z$distribution==
"power variance function Poisson")
z$coef[z$np] else exp(z$coef[z$np])))
colname <- c(colname,"parameter")}
dimnames(coef.table) <- list(cname,colname)}
print.default(coef.table,digits=digits,print.gap=2)}
if(z$archtype>0){
cat(if(z$archtype==1)"\nSquare"else if(z$archtype==2)"\nAbsolute value"
else "\nExponential","ARCH parameter\n")
cname <- c("estimate","se",if(z$archtype!=3)"arch")
coef.table <- cbind(z$coef[z$npa],z$se[z$npa],
if(z$archtype!=3)exp(z$coef[z$npa]))
dimnames(coef.table) <- list("", cname)
print.default(coef.table,digits=digits,print.gap=2)}
if(correlation){
cat("\nCorrelation matrix\n")
print.default(z$corr, digits=digits)}}
### volatility method
###
#' @export
volatility <- function(z, ...) UseMethod("volatility")
#' @export
volatility.gar <- function(z, nind=NULL, ...){
if(is.null(nind))nind <- 1:dim(z$response$y)[1]
else if(length(nind)>length(nobs(z))||any(nind>length(nobs(z))))
stop("Individual not found")
else nind <- !is.na(match(covind(z),nind))
if(all(!nind))stop("No such individuals")
if(is.null(z$volatility))return(NULL)
z$volatility[nind]}
Try the repeated package in your browser
Any scripts or data that you put into this service are public.
repeated documentation built on Aug. 8, 2023, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions volatility.gar volatility print.gar residuals.gar fitted.gar deviance.gar gar
Documented in deviance.gar fitted.gar gar print.gar residuals.gar volatility volatility.gar
#
# repeated : A Library of Repeated Measurements Models
# Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public Licence as published by
# the Free Software Foundation; either version 2 of the Licence, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public Licence for more details.
#
# You should have received a copy of the GNU General Public Licence
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
# SYNOPSIS
#
# gar(response=NULL, distribution="normal", times=NULL, totals=NULL,
# censor=NULL, delta=NULL, mu=NULL, shape=NULL, depend=NULL,
# shfn=FALSE, common=FALSE, preg=NULL, pshape=NULL, pdepend=NULL,
# parch=NULL, arch="square", transform="identity", link="identity",
# autocorr="exponential", order=1, envir=parent.frame(), print.level=0,
# ndigit=10, gradtol=0.00001, steptol=0.00001, fscale=1, iterlim=100,
# typsize=abs(p), stepmax=10*sqrt(p%*%p))
#
# DESCRIPTION
#
# A function to fit generalized nonlinear autoregression models with
# various distributions
##' Generalized Autoregression Models
##'
##' \code{gar} fits a first- or second-order generalized autoregression,
##' possibly with Kalman update over time (first-order only).
##'
##' Nonlinear regression models can be supplied as formulae where parameters
##' are unknowns in which case factor variables cannot be used and parameters
##' must be scalars. (See \code{\link[rmutil]{finterp}}.)
##'
##' Marginal and individual profiles can be plotted using
##' \code{\link[rmutil]{mprofile}} and \code{\link[rmutil]{iprofile}} and
##' residuals with \code{\link[rmutil]{plot.residuals}}.
##'
##' When the dispersion parameter is not constant over time, \code{volatility}
##' extracts the square root of the dispersion parameter for a fitted model.
##'
##'
##' @aliases gar volatility volatility.gar deviance.gar fitted.gar print.gar residuals.gar volatility.gar
##' @param response A list of two or three column matrices with responses,
##' corresponding times, and possibly a censor indicator, for each individual,
##' one matrix or dataframe of responses, or an object of class,
##' \code{response} (created by \code{\link[rmutil]{restovec}}) or
##' \code{repeated} (created by \code{\link[rmutil]{rmna}}). If the
##' \code{repeated} data object contains more than one response variable, give
##' that object in \code{envir} and give the name of the response variable to
##' be used here.
##' @param distribution The distribution to be fitted: binomial, Poisson,
##' exponential, negative binomial, mult Poisson, double Poisson, Consul
##' generalized Poisson, beta binomial, mult binomial, double binomial, normal,
##' inverse Gauss, logistic, gamma, Weibull, Cauchy, Laplace, Levy, Pareto,
##' beta, simplex, two-sided power, gen(eralized) gamma, gen(eralized)
##' logistic, Hjorth, Burr, gen(eralized) Weibull, gen(eralized) extreme value,
##' gen(eralized) inverse Gauss, power exponential, power variance function
##' Poisson, skew Laplace, or Student t. (For definitions of distributions, see
##' the corresponding [dpqr]distribution help.)
##' @param 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,
##' \code{response} or \code{repeated}.
##' @param totals An appropriate scalar, vector, or matrix of binomial totals
##' (only applicable for binomial, beta binomial, mult binomial, double
##' binomial). Ignored if response has class, \code{response} or
##' \code{repeated}.
##' @param censor If response is a matrix, a matrix of the same size containing
##' the censor indicator: 1=uncensored, 0=right-censored, -1=left-censored.
##' Ignored if response has class, \code{response} or \code{repeated}.
##' @param 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, \code{delta=0.01}. If the response has been
##' pretransformed, this must be multiplied by the Jacobian. This
##' transformation cannot contain unknown parameters. For example, with a log
##' transformation, \code{delta=1/y}. (The delta values for the censored
##' response are ignored.) The jacobian is calculated automatically for the
##' transform option. Ignored if response has class, \code{response} or
##' \code{repeated}.
##' @param mu A user-specified function of \code{pmu} giving the regression
##' equation for the location. It may also be 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. It must yield a value for each observation on each
##' individual.
##' @param shape An optional user-specified shape regression function; this may
##' depend on the location (function) through its second argument, in which
##' case, \code{shfn} must be set to TRUE. It may also be a formula beginning
##' with ~, specifying either a linear regression function for the shape
##' parameter in the Wilkinson and Rogers notation or a general function with
##' named unknown parameters. If it contains unknown parameters, the keyword
##' \code{mu} may be used to specify a function of the location parameter.
##' @param depend An optional user-specified regression function for the log
##' dependence parameter. It may also be a formula beginning with ~, specifying
##' either a linear regression function for the dependence parameter in the
##' Wilkinson and Rogers notation or a general function with named unknown
##' parameters. If used, \code{order} must be one.
##' @param shfn If TRUE, the supplied shape function depends on the location
##' function. The name of this location function must be the last argument of
##' the shape function.
##' @param common If TRUE, \code{mu} and \code{shape} must both be either
##' 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
##' or formulae with unknowns. All parameter estimates must be supplied in
##' \code{preg}. If FALSE, parameters are distinct between the two functions
##' and indexing starts at one in each function.
##' @param preg The initial parameter estimates for the location regression
##' function. If \code{mu} is a formula with unknown parameters, their
##' estimates must be supplied either in their order of appearance in the
##' expression or in a named list.
##' @param pdepend One or two estimates of the dependence parameters for the
##' Kalman update. With one, it is Markovian and, with two, it is
##' nonstationary. For the latter, the \code{order} must be one. If
##' \code{depend} is a function or formula, the corresponding number of
##' estimates must be supplied. Either pdepend or parch or both must be
##' supplied.
##' @param parch Estimate for an ARCH model where the shape parameter depends
##' on the square of the previous residual. Either pdepend or parch or both
##' must be supplied.
##' @param arch If \code{square}, then \code{shape+parch^diff*residual^2}; if
##' \code{absolute value}, then \code{shape+parch^diff*|residual|}; if
##' \code{exponential}, then \code{shape*exp(parch*residual^2*diff)}, where
##' \code{diff} is the length of time since the previous observation and
##' \code{residual} is the previous residual or innovation.
##' @param pshape Zero to two estimates for the shape parameters, depending on
##' the distribution, if \code{shape} is not a function; otherwise, estimates
##' for the parameters in this function, with one extra at the end for
##' three-parameter distributions. If \code{shape} is a formula with unknown
##' parameters, their estimates must be supplied either in their order of
##' appearance in the expression or in a named list (only for two-parameter
##' distributions).
##' @param transform Transformation of the response variable: \code{identity},
##' \code{exp}, \code{square}, \code{sqrt}, or \code{log}.
##' @param link Link function for the mean: \code{identity}, \code{exp},
##' \code{square}, \code{sqrt}, \code{log}, \code{logit}, \code{cloglog} or
##' \code{loglog} (last three only for binary/binomial-type data).
##' @param autocorr The form of the (second if two) dependence function:
##' \code{exponential} is the usual \eqn{\rho^{|t_i-t_j|}}{rho^|t_i-t_j|};
##' \code{gaussian} is \eqn{\rho^{(t_i-t_j)^2}}{rho^((t_i-t_j)^2)};
##' \code{cauchy} is \eqn{1/(1+\rho(t_i-t_j)^2)}{1/(1+rho(t_i-t_j)^2)};
##' \code{spherical} is
##' \eqn{((|t_i-t_j|\rho)^3-3|t_i-t_j|\rho+2)/2}{((|t_i-t_j|rho)^3-3|t_i-t_j|rho+2)/2}
##' for \eqn{|t_i-t_j|\leq1/\rho}{|t_i-t_j|<=1/rho} and zero otherwise;
##' \code{IOU} is the integrated Ornstein-Uhlenbeck process, \eqn{(2\rho
##' \min(t_i,t_j)+\exp(-\rho t_i) }{(2rho min(t_i,t_j)+exp(-rho t_i) +exp(-rho
##' t_j)-1 -exp(rho|ti-t_j|))/2rho^3}\eqn{+\exp(-\rho t_j)-1
##' -\exp(\rho|ti-t_j|))/2\rho^3}{(2rho min(t_i,t_j)+exp(-rho t_i) +exp(-rho
##' t_j)-1 -exp(rho|ti-t_j|))/2rho^3}.
##' @param order First- or second-order stationary autoregression.
##' @param envir Environment in which model formulae are to be interpreted or a
##' data object of class, \code{repeated}, \code{tccov}, or \code{tvcov}; the
##' name of the response variable should be given in \code{response}. If
##' \code{response} has class \code{repeated}, it is used as the environment.
##' @param print.level Arguments for nlm.
##' @param typsize Arguments for nlm.
##' @param ndigit Arguments for nlm.
##' @param gradtol Arguments for nlm.
##' @param stepmax Arguments for nlm.
##' @param steptol Arguments for nlm.
##' @param iterlim Arguments for nlm.
##' @param fscale Arguments for nlm.
##' @return A list of classes \code{gar} and \code{recursive} is returned that
##' contains all of the relevant information calculated, including error codes.
##'
##' The volatility vector for models with a shape regression function and ARCH
##' models contains the square root of the dispersion parameter at each time
##' point.
##' @author J.K. Lindsey
### @seealso \code{\link[growth]{carma}}, \code{\link[growth]{elliptic}},
### \code{\link[rmutil]{finterp}}, \code{\link[repeated]{gnlmm}},
### \code{\link[gnlm]{gnlr}}, \code{\link[rmutil]{iprofile}},
### \code{\link[repeated]{kalcount}}, \code{\link[repeated]{kalseries}},
### \code{\link[event]{kalsurv}}, \code{\link[rmutil]{mprofile}},
### \code{\link[rmutil]{plot.residuals}}, \code{\link[rmutil]{read.list}},
### \code{\link[rmutil]{restovec}}, \code{\link[rmutil]{rmna}},
### \code{\link[rmutil]{tcctomat}}, \code{\link[rmutil]{tvctomat}}.
##' @references Lindsey, J.K. (1997) Applying Generalized Linear Models.
##' Springer, pp.\ 93--101
##'
##' Lambert, P. (1996) Statistics in Medicine 15, 1695-1708
##' @keywords models
##' @examples
##'
##' # first-order one-compartment model
##' # data objects for formulae
##' dose <- c(2,5)
##' dd <- tcctomat(dose)
##' times <- matrix(rep(1:20,2), nrow=2, byrow=TRUE)
##' tt <- tvctomat(times)
##' # vector covariates for functions
##' dose <- c(rep(2,20),rep(5,20))
##' times <- rep(1:20,2)
##' # functions
##' 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)
##' # response
##' conc <- matrix(rgamma(40,shape(log(c(0.1,0.4))),
##' scale=mu(log(c(1,0.3,0.2))))/shape(log(c(0.1,0.4))),ncol=20,byrow=TRUE)
##' conc[,2:20] <- conc[,2:20]+0.5*(conc[,1:19]-matrix(mu(log(c(1,0.3,0.2))),
##' ncol=20,byrow=TRUE)[,1:19])
##' conc <- restovec(ifelse(conc>0,conc,0.01),name="conc")
##' reps <- rmna(conc, ccov=dd, tvcov=tt)
##' # constant shape parameter
##' gar(conc, dist="gamma", times=1:20, mu=mu,
##' preg=log(c(1,0.4,0.1)), pdepend=0.5, pshape=1)
##' \dontrun{# or
##' gar(conc, dist="gamma", times=1:20, mu=~exp(absorption-volume)*
##' dose/(exp(absorption)-exp(elimination))*
##' (exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
##' preg=list(absorption=1,elimination=log(0.4),volume=log(0.1)),
##' pdepend=0.5, pshape=1, envir=reps)
##' # generalized gamma distribution
##' gar(conc, dist="gen gamma", times=1:20, mu=mu,
##' preg=log(c(1,0.4,0.1)), pdepend=0.3, pshape=c(.1,1))
##' # (if the covariates contained NAs, reps would have to be used as
##' # response instead of conc)
##' #
##' # time dependent shape parameter
##' gar(conc, dist="gamma", times=1:20, mu=mu, shape=shape,
##' preg=log(c(1,0.4,0.1)), pdepend=0.25, pshape=c(exp(-2),exp(-.57)))
##' # or
##' gar(conc, dist="gamma", times=1:20, mu=~exp(absorption-volume)*
##' dose/(exp(absorption)-exp(elimination))*
##' (exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
##' shape=~exp(b1-b2)*times*dose*exp(-exp(b1)*times),
##' preg=list(absorption=0,elimination=log(0.4),volume=log(0.1)),
##' pdepend=0.3, pshape=list(b1=exp(-2),b2=exp(-.57)), envir=reps)
##' # generalized gamma distribution
##' gar(conc, dist="gen gamma", times=1:20, mu=mu, shape=shape,
##' preg=log(c(1,0.4,0.1)), pdepend=0.5,
##' pshape=c(exp(-2),exp(-.57),2))
##' #
##' # shape function depends on location parameter
##' shape <- function(p, mu) p[1]+p[2]*mu
##' gar(conc, dist="gamma", times=1:20, mu=mu, shape=shape, shfn=TRUE,
##' preg=log(c(1,0.4,.10)), pdepend=0.15, pshape=c(1,2))
##' # or
##' gar(conc, dist="gamma", times=1:20, mu=mu, shape=~a+d*mu, shfn=TRUE,
##' preg=log(c(1,0.4,.10)), pdepend=0.15, pshape=c(1,2))
##'
##' }
##' @export gar
gar <- function(response=NULL, distribution="normal", times=NULL, totals=NULL,
censor=NULL, delta=NULL, mu=NULL, shape=NULL, depend=NULL, shfn=FALSE,
common=FALSE, preg=NULL, pshape=NULL, pdepend=NULL, parch=NULL,
arch="square", transform="identity", link="identity",
autocorr="exponential", order=1, envir=parent.frame(), print.level=0,
ndigit=10, gradtol=0.00001, steptol=0.00001, fscale=1, iterlim=100,
typsize=abs(p), stepmax=10*sqrt(p%*%p)){
#
# set up likelihood to call C code
#
likekal <- function(p){
eta <- mu1(p)
if(sh){
shr <- sh1(p)
theta <- exp(p[np])
if(!tm)theta <- c(p[npr1:nprd],theta)}
else theta <- if(tm)p[nprd1:np] else p[npr1:np]
if(tm)dep <- dep1(p)
if(archtype>0)parch <- p[npa]
z <- .C("gar_c",
y=as.double(y),
total=as.double(response$response$n),
my=as.integer(3*max(y)),
nobs=as.integer(nobs(response)),
nind=as.integer(nind),
times=as.double(response$response$times),
censor=as.integer(censor),
cens=as.integer(!is.null(response$response$censor)),
eta=as.double(eta),
tm=as.integer(tm),
theta=as.double(theta),
dep=as.double(dep),
model=as.integer(mdl),
thp=as.integer(thp),
shape=as.double(shr),
sh=as.integer(sh),
link=as.integer(lnk),
ar=as.integer(ar),
order=as.integer(order),
parch=as.double(parch),
arch=as.integer(archtype),
pred=double(n),
rpred=double(n),
volatility=double(n),
like=double(1),
#DUP=FALSE,
PACKAGE="repeated")
z$like+jacob}
call <- sys.call()
#
# check model required
#
tmp <- c("binomial","Poisson","exponential","negative binomial",
"mult Poisson","double Poisson","Consul","beta binomial",
"mult binomial","double binomial","normal","logistic","Cauchy",
"Weibull","gamma","Laplace","inverse Gauss","Pareto","Levy",
"beta","simplex","two-sided power","gen gamma","gen logistic","Hjorth",
"Burr","gen Weibull","gen extreme value","gen inverse Gauss",
"power exponential","power variance function Poisson",
"skew Laplace","Student t")
mdl <- match(distribution <- match.arg(distribution,tmp),tmp)
bmdl <- mdl==1||mdl==8||mdl==9||mdl==10
par1 <- distribution=="binomial"||distribution=="Poisson"||
distribution=="exponential"
par3 <- mdl>22
tmp <- c("exponential","gaussian","cauchy","spherical","IOU")
ar <- match(autocorr <- match.arg(autocorr,tmp),tmp)
#
# check transformation and link function
#
transform <- match.arg(transform,c("identity","exp","square","sqrt","log"))
if(transform!="identity"&&((mdl<11&&mdl!=3)||mdl==31))
stop("transformations not valid for discrete data models")
tmp <- c("identity","exp","square","sqrt","log","logit","cloglog","loglog")
lnk <- match(link <- match.arg(link,tmp),tmp)
if((link=="logit"||link=="cloglog"||link=="loglog")&&!bmdl)
stop("logit, cloglog, and loglog links can only be used with binary data")
#
# check that either AR or ARCH is being fitted
#
if(is.null(pdepend)&&is.null(parch))
stop("Either pdepend or parch must be supplied")
if(!is.null(parch)){
if(par1)stop(paste("arch model not possible for",distribution,"distribution"))
if(length(parch)!=1)stop("parch should contain one value")
tmp <- c("square","absolute value","exponential")
archtype <- match(match.arg(arch,tmp),tmp)
if(archtype!=3&&parch<=0)
stop("parch must be positive for square and absolute value models")}
else archtype <- 0
#
# check if a data object is being supplied
#
type <- "unknown"
respenv <- exists(deparse(substitute(response)),envir=parent.frame())&&
inherits(response,"repeated")
if(respenv){
if(dim(response$response$y)[2]>1)
stop("gar only handles univariate responses")
if(!is.null(response$NAs)&&any(response$NAs))
stop("gar does not handle data with NAs")
type <- response$response$type}
envname <- if(respenv)deparse(substitute(response))
else if(!is.null(class(envir)))deparse(substitute(envir))
else NULL
#
#if response is not a data object, find or make one
#
if(!respenv){
if(inherits(envir,"repeated")){
# if envir, remove extra (multivariate) responses
if(!is.null(envir$NAs)&&any(envir$NAs))
stop("gar does not handle data with NAs")
cn <- deparse(substitute(response))
if(length(grep("\"",cn))>0)cn <- response
if(length(cn)>1)stop("only one response variable allowed")
response <- envir
col <- match(cn,colnames(response$response$y))
if(is.na(col))stop(paste("response variable",cn,"not found"))
type <- response$response$type[col]
if(dim(response$response$y)[2]>1){
response$response$y <- response$response$y[,col,drop=FALSE]
if(!is.null(response$response$delta)){
response$response$delta <- response$response$delta[,col,drop=FALSE]
if(all(response$response$delta==1)||all(is.na(response$response$delta)))response$response$delta <- NULL}
if(!is.null(response$response$n)){
response$response$n <- response$response$n[,col,drop=FALSE]
if(all(is.na(response$response$n)))response$response$n <- NULL}
if(!is.null(envir$response$censor)){
response$response$censor <- response$response$censor[,col,drop=FALSE]
if(all(is.na(response$response$censor)))response$response$censor <- NULL}}}
else {
if(!inherits(response,"response")){
if(is.null(times)){
if(is.matrix(response))times <- matrix(1:dim(response)[2],nrow=dim(response)[1],ncol=dim(response)[2],byrow=T)
else if(is.vector(response))times <- 1:length(response)}
response <- if(bmdl)
restovec(response, times=times, totals=totals)
else restovec(response, times=times, censor=censor, delta=delta)}
type <- response$type
response <- rmna(response=response)}}
if(is.null(times(response)))stop("No times available")
if((inherits(envir,"repeated")&&
(length(nobs(response))!=length(nobs(envir))||
any(nobs(response)!=nobs(envir))))||(inherits(envir,"tvcov")&&
(length(nobs(response))!=length(envir$tvcov$nobs)||
any(nobs(response)!=envir$tvcov$nobs))))
stop("response and envir objects are incompatible")
if(bmdl){
if(type!="unknown"&&type!="nominal")stop("nominal data required")
if(is.null(response$response$n)){
if(any(response$response$y!=0&&response$response$y!=1))
stop("responses must be binary if totals are not supplied")
else response$response$n <- rep(1,dim(response$response$y)[1])}}
y <- response$response$y
n <- dim(y)[1]
#
# if a data object was supplied, modify formulae or functions to read from it
#
mu3 <- sh3 <- dep3 <- dep <- NULL
if(respenv||inherits(envir,"repeated")||inherits(envir,"tccov")||inherits(envir,"tvcov")){
if(is.null(envname))envname <- deparse(substitute(envir))
# location function or formula
if(inherits(mu,"formula")){
mu3 <- if(respenv)finterp(mu,.envir=response,.name=envname)
else finterp(mu,.envir=envir,.name=envname)}
else if(is.function(mu)){
if(is.null(attr(mu,"model"))){
tmp <- parse(text=deparse(mu)[-1])
mu <- if(respenv)fnenvir(mu,.envir=response,.name=envname)
else fnenvir(mu,.envir=envir,.name=envname)
mu3 <- mu
attr(mu3,"model") <- tmp}
else mu3 <- mu}
# shape function or formula
if(inherits(shape,"formula")){
sh3 <- if(respenv)finterp(shape,.envir=response,.name=envname)
else finterp(shape,.envir=envir,.name=envname)}
else if(is.function(shape)){
if(is.null(attr(shape,"model"))){
tmp <- parse(text=deparse(shape)[-1])
shape <- if(respenv)fnenvir(shape,.envir=response,.name=envname)
else fnenvir(shape,.envir=envir,.name=envname)
sh3 <- shape
if(shfn){
pos <- grep("mu",attr(sh3,"parameters"))
attr(sh3,"parameters") <- attr(sh3,"parameters")[-pos]}
attr(sh3,"model") <- tmp}
else sh3 <- shape}
# AR function or formula
if(inherits(depend,"formula")){
dep3 <- if(respenv)finterp(depend,.envir=response,.name=envname)
else finterp(depend,.envir=envir,.name=envname)}
else if(is.function(depend)){
if(is.null(attr(depend,"model"))){
tmp <- parse(text=deparse(depend)[-1])
depend <- if(respenv)fnenvir(depend,.envir=response,.name=envname)
else fnenvir(depend,.envir=envir,.name=envname)
dep3 <- depend
attr(dep3,"model") <- tmp}
else dep3 <- depend}}
#
# transform location formula to function and check number of parameters
#
npr <- length(preg)
npr1a <- npr1 <- npr+1
if(inherits(mu,"formula")){
mu2 <- if(respenv)finterp(mu,.envir=response,.name=envname)
else finterp(mu,.envir=envir,.name=envname)
npt1 <- length(attr(mu2,"parameters"))
if(is.character(attr(mu2,"model"))){
# W&R formula
if(length(attr(mu2,"model"))==1){
# intercept model
mu1 <- function(p) p[1]*rep(1,n)
attributes(mu1) <- attributes(mu2)
mu2 <- NULL}}
else {
# formula with unknowns
if(npr!=npt1&&!common){
cat("\nParameters are ")
cat(attr(mu2,"parameters"),"\n")
stop(paste("preg should have",npt1,"estimates"))}
if(is.list(preg)){
if(!is.null(names(preg))){
o <- match(attr(mu2,"parameters"),names(preg))
preg <- unlist(preg)[o]
if(sum(!is.na(o))!=length(preg))stop("invalid estimates for mu - probably wrong names")}
else preg <- unlist(preg)}}
if(!is.null(mu2)){
if(inherits(envir,"tccov")){
cv <- covind(response)
mu1 <- function(p) mu2(p)[cv]
attributes(mu1) <- attributes(mu2)}
else {
mu1 <- mu2
rm(mu2)}}}
else if(is.null(mu)){
mu1 <- function(p) p[1]*rep(1,n)
npt1 <- 1}
else mu1 <- mu
#
# give appropriate attributes to mu1 for printing
#
if(!is.null(mu1)&&is.null(attr(mu1,"parameters"))){
attributes(mu1) <- if(is.function(mu)){
if(!inherits(mu,"formulafn")){
if(respenv)attributes(fnenvir(mu,.envir=response))
else attributes(fnenvir(mu,.envir=envir))}
else attributes(mu)}
else {
if(respenv)attributes(fnenvir(mu1,.envir=response))
else attributes(fnenvir(mu1,.envir=envir))}}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(mu1))length(attr(mu1,"parameters"))
else if(is.null(mu1))NULL
else npt1
if(!is.null(nlp)&&!common&&nlp!=npr)
stop(paste("preg should have",nlp,"initial estimates"))
#
# transform shape formula to function and check number of parameters
#
nprd <- npr+length(pdepend)
nprd1 <- if(common&&!inherits(shape,"formula")) 1 else nprd+1
nps1 <- nps <- length(pshape)
if(par3)nps1 <- nps1-1
if(inherits(shape,"formula")){
old <- if(common)mu1 else NULL
mufn <- if(shfn)"mu" else NULL
sh3 <- if(respenv)finterp(shape,.envir=response,.start=nprd1,.name=envname,.old=old,.args=mufn)
else finterp(shape,.envir=envir,.start=nprd1,.name=envname,.old=old,.args=mufn)
tmp <- attributes(sh3)
sh2 <- if(shfn)function(p) sh3(p,mu1(p)) else sh3
attributes(sh2) <- tmp
npt2 <- length(attr(sh2,"parameters"))
if(is.character(attr(sh2,"model"))){
# W&R formula
if(length(attr(sh2,"model"))==1){
# intercept model
sh1 <- function(p) p[nprd1]*rep(1,n)
attributes(sh1) <- attributes(sh2)
sh2 <- NULL}}
else {
# formula with unknowns
if(nps1!=npt2&&!common){
cat("\nParameters are ")
cat(attr(sh2,"parameters"),"\n")
stop(paste("pshape should have",npt2+par3,"estimates"))}
if(is.list(pshape)){
if(!is.null(names(pshape))){
o <- match(attr(sh2,"parameters"),names(pshape))
pshape <- unlist(pshape)[o]
if(sum(!is.na(o))!=length(pshape))stop("invalid estimates for shape - probably wrong names")}
else pshape <- unlist(pshape)}}
if(!is.null(sh2)){
if(inherits(envir,"tccov")){
cv <- covind(response)
sh1 <- function(p) sh2(p)[cv]
attributes(sh1) <- attributes(sh2)}
else {
sh1 <- sh2
rm(sh2)}}}
else if(!is.function(shape)&&!par1){
sh1 <- function(p) p[nprd1]*rep(1,n)
npt2 <- 1
if(length(pshape)!=1+par3)
stop(paste("pshape must contain",1+par3,"estimate(s)"))}
else sh1 <- if(shfn)function(p) shape(p[nprd1:np], mu1(p))
else function(p) shape(p[nprd1:np])
#
# give appropriate attributes to sh1 for printing
#
if(!is.null(sh1)&&is.null(attr(sh1,"parameters")))
attributes(sh1) <- if(is.function(shape)){
if(!inherits(shape,"formulafn")){
if(respenv)attributes(fnenvir(shape,.envir=response))
else attributes(fnenvir(shape,.envir=envir))}
else attributes(shape)}
else {
if(respenv)attributes(fnenvir(sh1,.envir=response))
else attributes(fnenvir(sh1,.envir=envir))}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(shape))length(attr(sh1,"parameters"))-shfn
else if(is.null(shape))NULL
else npt2
if(!is.null(nlp)&&!common&&nlp!=nps1)
stop(paste("pshape should have",nlp+par3,"initial estimates"))
if(common){
nlp <- length(unique(c(attr(mu1,"parameters"),attr(sh1,"parameters"))))
if(nlp!=npr)stop(paste("with a common parameter model, preg should contain",nlp,"estimates"))}
#
# transform AR formula to function and check number of parameters
#
npd <- length(pdepend)
if(inherits(depend,"formula")){
dep2 <- if(respenv)finterp(depend,.envir=response,.start=npr1,.name=envname)
else finterp(depend,.envir=envir,.start=npr1,.name=envname)
npt3 <- length(attr(dep2,"parameters"))
if(is.character(attr(dep2,"model"))){
# W&R formula
if(length(attr(dep2,"model"))==1){
# intercept model
dep1 <- function(p) p[npr1]*rep(1,n)
attributes(dep1) <- attributes(dep2)
dep2 <- NULL}}
else {
# formula with unknowns
if(npd!=npt3){
cat("\nParameters are ")
cat(attr(dep2,"parameters"),"\n")
stop(paste("pdepend should have",npt3,"estimates"))}
if(is.list(pdepend)){
if(!is.null(names(pdepend))){
o <- match(attr(dep2,"parameters"),names(pdepend))
pdepend <- unlist(pdepend)[o]
if(sum(!is.na(o))!=length(pdepend))stop("invalid estimates for depend - probably wrong names")}
else pdepend <- unlist(pdepend)}}
if(!is.null(dep2)){
if(inherits(envir,"tccov")){
cv <- covind(response)
dep1 <- function(p) dep2(p)[cv]
attributes(dep1) <- attributes(dep2)}
else {
dep1 <- dep2
rm(dep2)}}}
else if(is.function(dep))dep1 <- dep
else {
dep1 <- NULL
npt3 <- 1}
#
# give appropriate attributes to dep1 for printing
#
if(!is.null(dep1)&&is.null(attr(dep1,"parameters")))
attributes(dep1) <- if(is.function(depend)){
if(!inherits(depend,"formulafn")){
if(respenv)attributes(fnenvir(depend,.envir=response))
else attributes(fnenvir(depend,.envir=envir))}
else attributes(depend)}
else {
if(respenv)attributes(fnenvir(dep1,.envir=response))
else attributes(fnenvir(dep1,.envir=envir))}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(depend))length(attr(dep1,"parameters"))
else if(is.null(depend))NULL
else npt3
if(!is.null(nlp)&&nlp!=npd)
stop(paste("pdepend should have",nlp,"initial estimates"))
#
# check response and functions
#
if(mdl<11){
if((distribution=="Poisson"||distribution=="negative binomial"
||distribution=="mult Poisson"
||distribution=="double Poisson"||distribution=="Consul")
&&type!="unknown"&&type!="discrete")
stop("discrete data required")
if(distribution=="exponential"&&type!="unknown"&&type!="duration")
stop("duration data required")
if(any(y<0))stop("All responses must be non-negative")}
else if(mdl==20||mdl==21||mdl==22){
if(type!="unknown"&&type!="continuous")stop("continuous data required")
if(any(y<=0)||any(y>=1))stop("All responses must be between 0 and 1")}
else if(mdl==31){
if(type!="unknown"&&type!="discrete")stop("discrete data required")
if(any(y<0))stop("All responses must be non-negative")}
else if(mdl!=11&&mdl!=12&&mdl!=13&&mdl!=16&&mdl!=19&&mdl!=24&&mdl!=30&&
mdl!=32&&mdl!=33){
if(type!="unknown"&&type!="duration"&&type!="continuous")
stop("duration data required")
if(any(y<=0))stop("All responses must be positive")}
else {
if(type!="unknown"&&type!="continuous"&&type!="duration")
stop("continuous data required")
if(distribution=="Levy"&&any(response$response$y<=mu1(preg)))
stop("Location function must give values strictly less than corresponding observations")}
if(distribution=="Pareto"&&any(sh1(pshape)<=1))
stop("shape function must be > 1")
#
# set up censoring indicator
#
censor <- response$response$censor
if(is.null(censor))censor <- rep(1,n)
else if(bmdl||mdl==2||mdl==4||mdl==5||mdl==6||mdl==7||mdl==31)
stop(paste("Censored data not allowed for",distribution,"distribution"))
#
# set up transformation and Jacobian
#
nind <- length(nobs(response))
if(transform=="identity")jacob <- 0
else if(transform=="exp"){
jacob <- -sum(y[censor==1])
y <- exp(y)}
else if(any(y==0))stop("Zero response values: invalid transformation")
else if(transform=="square"){
jacob <- -sum(log(abs(y[censor==1])))-n*log(2)
y <- y^2}
else if(any(y<0))stop("Nonpositive response values: invalid transformation")
else if(transform=="sqrt"){
jacob <- sum(log(y[censor==1]))/2+n*log(2)
y <- sqrt(y)}
else if(transform=="log"){
jacob <- sum(log(y[censor==1]))
y <- log(y)}
#
# set up unit of measurement
#
if(!is.null(response$response$delta)){
if(length(response$response$delta)==1)
jacob <- jacob-length(y[censor==1])*log(response$response$delta)
else jacob <- jacob-sum(log(response$response$delta[censor==1]))}
#
# check AR model
#
if(is.null(pdepend))order <- 0
else if(order!=1&&order!=2)stop("Autoregression must have order 1 or 2")
tm <- !is.null(depend)
if(tm)order <- 0
if(order==2&&length(pdepend)!=2)
stop("2 estimates of dependence parameters must be given")
else if(order>0&&!tm&&length(pdepend)!=1&&length(pdepend)!=2)
stop("One or two estimates of dependence parameters must be given")
thp <- length(pdepend)==2&&order==1&&!tm
if(!tm&&any(pdepend<=0))stop("All dependence parameters must be positive")
#
# check location function
#
if(length(mu1(preg))!=dim(response$response$y)[1])
stop("The mu function must provide an estimate for each observation")
else if(any(is.na(mu1(preg))))
stop("Non-numerical mu: probably invalid initial values")
#
# set up parameter vector
#
p <- if(tm)c(preg,pdepend) else if(order>0)c(preg,-log(pdepend)) else preg
sh <- mdl>3&&(is.function(shape)||inherits(shape,"formula"))
if(!par1){
if(!sh)p <- c(p,if(mdl==31)c(log(pshape[1]),pshape[2])else log(pshape))
else {
if(par3){
if(pshape[nps]<=0&&mdl!=31)
stop("extra shape parameter must be positive")
p <- if(common)c(p,if(mdl==31)pshape[nps]
else log(pshape[nps]))
else c(p,pshape[1:(nps-1)],if(mdl==31)
pshape[nps] else log(pshape[nps]))}
else p <- c(p,pshape)}}
np <- length(p)
if(archtype>0)p <- c(p,if(archtype==1)log(parch)else parch)
npa <- length(p)
#
# check shape function
#
if(!sh){
if((mdl<=3&&nps!=0)||(mdl>3&&mdl<23&&nps!=1)||(mdl>=23&&nps!=2))
stop("Incorrect number of shape parameter estimates")
else if(nps>0){
if(mdl!=31&&any(pshape<=0)||mdl==31&&pshape[1]<=0)
stop("All shape parameters must be positive")}
shr <- rep(0,dim(response$response$y)[1])}
else {
if(any(is.na(sh1(p))))
stop("The shape model returns NAs: probably invalid initial values")
if(length(sh1(p))!=dim(response$response$y)[1])
stop("The shape function must provide an estimate for each observation")}
#
# check likelihood function and optimize with nlm
#
if(fscale==1)fscale <- likekal(p)
if(is.na(likekal(p)))
stop("Likelihood returns NAs: probably invalid initial values")
z0 <- nlm(likekal, p, hessian=TRUE, print.level=print.level,
typsize=typsize, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
steptol=steptol, iterlim=iterlim, fscale=fscale)
like <- z0$minimum
#
# calculate se's
#
a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
else qr(z0$hessian)$rank
if(a==npa)cov <- solve(z0$hessian)
else cov <- matrix(NA,ncol=npa,nrow=npa)
se <- sqrt(diag(cov))
corr <- cov/(se%o%se)
dimnames(corr) <- list(1:npa,1:npa)
#
# back transform parameters
#
eta <- mu1(z0$estimate)
if(sh){
shr <- sh1(z0$estimate)
theta <- exp(z0$estimate[np])
if(!tm)theta <- c(z0$estimate[npr1:nprd],theta)}
else theta <- if(tm)z0$estimate[nprd1:np] else z0$estimate[npr1:np]
if(tm)dep <- dep1(z0$estimate)
#
# calculate predicted values (and volatility)
#
if(archtype>0)parch <- z0$estmate[npa]
z <- .C("gar_c",
y=as.double(y),
total=as.double(response$response$n),
my=as.integer(3*max(y)),
nobs=as.integer(nobs(response)),
nind=as.integer(nind),
times=as.double(response$response$times),
censor=as.integer(censor),
cens=as.integer(!is.null(response$response$censor)),
eta=as.double(eta),
tm=as.integer(tm),
theta=as.double(theta),
dep=as.double(dep),
model=as.integer(mdl),
thp=as.integer(thp),
shape=as.double(shr),
sh=as.integer(sh),
link=as.integer(lnk),
ar=as.integer(ar),
order=as.integer(order),
parch=as.double(parch),
arch=as.integer(archtype),
pred=double(n),
rpred=double(n),
volatility=double(n),
like=double(1),
#DUP=FALSE,
PACKAGE="repeated")
if(archtype==0)z$volatility <- NULL
#
# return object
#
if(!is.null(mu3))mu1 <- mu3
if(!is.null(sh3))sh1 <- sh3
if(!is.null(dep3))dep1 <- dep3
z <- list(
call=call,
distribution=distribution,
mu=mu1,
formula=mu,
shape=shape,
depend=depend,
sh1=sh1,
shfn=shfn,
common=common,
dep1=dep1,
response=response$response,
link=link,
order=order,
autocorr=autocorr,
transform=transform,
maxlike=like,
aic=like+npa,
df=dim(response$response$y)[1]-npa,
np=np,
npa=npa,
npr=npr,
nps=nps,
npd=npd,
thp=thp,
tm=tm,
archtype=archtype,
coefficients=z0$estimate,
se=se,
cov=cov,
corr=corr,
pred=z$pred,
rpred=z$rpred,
volatility=z$volatility,
grad=z0$gradient,
iterations=z0$iterations,
code=z0$code)
class(z) <- c("gar","recursive")
return(z)}
### standard methods
###
#' @export
deviance.gar <- function(object, ...) 2*object$maxlike
#' @export
fitted.gar <- function(object, recursive=TRUE, ...) if(recursive) object$rpred else object$pred
#' @export
residuals.gar <- function(object, recursive=TRUE, ...){
z <- object
if(z$transform=="exp")z$response$y <- exp(z$response$y)
else if(z$transform=="square")z$response$y <- z$response$y^2
else if(z$transform=="sqrt")z$response$y <- sqrt(z$response$y)
else if(z$transform=="log")z$response$y <- log(z$response$y)
if(recursive) z$response$y-z$rpred else z$response$y-z$pred}
### print method
###
#' @export
print.gar <- function(x,digits=max(3,.Options$digits-3),correlation=TRUE, ...){
z <- x
np1 <- if(z$distribution=="binomial"||z$distribution=="exponential"
||z$distribution=="Poisson") 0
else if(z$distribution=="gen gamma"
||z$distribution=="gen logistic"
||z$distribution=="Hjorth"||z$distribution=="Burr"
||z$distribution=="gen Weibull"
||z$distribution=="gen extreme value"
||z$distribution=="gen inverse Gauss"
||z$distribution=="power exponential"
||z$distribution=="power variance function Poisson"
||z$distribution=="skew Laplace"
||z$distribution=="Student t") 2
else 1
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
cat("Number of subjects ",length(nobs(z)),"\n")
cat("Number of observations",length(z$response$y),"\n")
cat("Transformation ",z$trans,"\n")
cat("Link function ",z$link,"\n\n")
cat(z$distribution,"distribution\n")
if(z$order<2)cat("First order ")
else cat("Second order ")
cat(z$autocorr,"dependence\n")
cat("\n-Log likelihood ",z$maxlike,"\n")
cat("Degrees of freedom",z$df,"\n")
cat("AIC ",z$aic,"\n")
cat("Iterations ",z$iterations,"\n\n")
if(z$common)cat("Location model\n")
else cat("Location parameters\n")
if(!is.null(attr(z$mu,"formula")))
cat(deparse(attr(z$mu,"formula")),sep="\n")
else if(!is.null(attr(z$mu,"model"))){
t <- deparse(attr(z$mu,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- if(is.character(attr(z$mu,"model")))attr(z$mu,"model")
else attr(z$mu,"parameters")
coef.table <- cbind(z$coef[1:z$npr],z$se[1:z$npr])
colname <- c("estimate","se")
if(!z$common){
dimnames(coef.table) <- list(cname,colname)
print.default(coef.table, digits=digits, print.gap=2)}
else {
cat("\nShape model\n")
if(!is.null(attr(z$sh1,"formula")))
cat(deparse(attr(z$sh1,"formula")),sep="\n")
else if(!is.null(attr(z$sh1,"model"))){
t <- deparse(attr(z$sh1,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- c(cname,if(is.character(attr(z$sh1,"model")))
attr(z$sh1,"model")
else attr(z$sh1,"parameters")[1:(length(attr(z$sh1,"parameters"))-z$shfn)])
cname <- unique(cname)
cat("\nCommon parameters\n")
dimnames(coef.table) <- list(cname,colname)
print.default(coef.table, digits=digits, print.gap=2)}
if(z$order>0||z$tm){
if(z$thp||z$order==2){
cat("\nDependence parameters\n")
if(z$thp)cname <- c("phi","rho")
else cname <- c("rho1","rho2")
coef.table <- cbind(z$coef[(z$npr+1):(z$npr+2)],
z$se[(z$npr+1):(z$npr+2)],
exp(-z$coef[(z$npr+1):(z$npr+2)]))}
else if(z$tm){
cat("\nDependence parameters\n")
if(!is.null(attr(z$dep1,"formula")))
cat(deparse(attr(z$dep1,"formula")),sep="\n")
else if(!is.null(attr(z$dep1,"model"))){
t <- deparse(attr(z$dep1,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- if(is.character(attr(z$dep1,"model")))
attr(z$dep1,"model")
else attr(z$dep1,"parameters")
coef.table <- cbind(z$coef[(z$npr+1):(z$npr+z$npd)],
z$se[(z$npr+1):(z$npr+z$npd)])
colname <- c("estimate","se")
dimnames(coef.table) <- list(cname,colname)}
else {
cat("\nDependence parameter\n")
cname <- "rho"
coef.table <- cbind(z$coef[z$npr+1],
z$se[z$npr+1],exp(-z$coef[z$npr+1]))}
if(!z$tm)dimnames(coef.table) <- list(cname, c("estimate","se","parameter"))
print.default(coef.table, digits=digits, print.gap=2)}
if(np1>0&&!z$common){
cat("\nShape parameters\n")
if(is.null(z$shape)){
cname <- "shape"
if(np1==2)cname <- c(cname,"psi")
coef.table <- cbind(z$coef[(z$np-np1+1):z$np],
z$se[(z$np-np1+1):z$np],if(z$distribution==
"power variance function Poisson")
c(exp(z$coef[(z$np-np1+1):(z$np-1)]),z$coef[z$np]) else
exp(z$coef[(z$np-np1+1):z$np]))
dimnames(coef.table) <- list(cname, c("estimate","se","parameter"))}
else {
if(!is.null(attr(z$sh1,"formula")))
cat(deparse(attr(z$sh1,"formula")),sep="\n")
else if(!is.null(attr(z$sh1,"model"))){
t <- deparse(attr(z$sh1,"model"))
t[1] <- sub("expression\\(","",t[1])
t[length(t)] <- sub("\\)$","",t[length(t)])
cat(t,sep="\n")}
cname <- if(is.character(attr(z$sh1,"model")))
attr(z$sh1,"model")
else attr(z$sh1,"parameters")
np2 <- length(cname)+(np1==2)
coef.table <- cbind(z$coef[(z$np-np2+1):z$np],
z$se[(z$np-np2+1):z$np])
if(np1==2){
cname <- c(cname,"psi")
coef.table <- cbind(coef.table,
c(coef.table[1:(dim(coef.table)[1]-1),1],
if(z$distribution==
"power variance function Poisson")
z$coef[z$np] else exp(z$coef[z$np])))
colname <- c(colname,"parameter")}
dimnames(coef.table) <- list(cname,colname)}
print.default(coef.table,digits=digits,print.gap=2)}
if(z$archtype>0){
cat(if(z$archtype==1)"\nSquare"else if(z$archtype==2)"\nAbsolute value"
else "\nExponential","ARCH parameter\n")
cname <- c("estimate","se",if(z$archtype!=3)"arch")
coef.table <- cbind(z$coef[z$npa],z$se[z$npa],
if(z$archtype!=3)exp(z$coef[z$npa]))
dimnames(coef.table) <- list("", cname)
print.default(coef.table,digits=digits,print.gap=2)}
if(correlation){
cat("\nCorrelation matrix\n")
print.default(z$corr, digits=digits)}}
### volatility method
###
#' @export
volatility <- function(z, ...) UseMethod("volatility")
#' @export
volatility.gar <- function(z, nind=NULL, ...){
if(is.null(nind))nind <- 1:dim(z$response$y)[1]
else if(length(nind)>length(nobs(z))||any(nind>length(nobs(z))))
stop("Individual not found")
else nind <- !is.na(match(covind(z),nind))
if(all(!nind))stop("No such individuals")
if(is.null(z$volatility))return(NULL)
z$volatility[nind]}
Try the repeated package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.