treat <- c(0,0,1,1)
tr <- tcctomat(treat)
dose <- matrix(rpois(20,10), ncol=5)
dd <- tvctomat(dose)
y <- restovec(matrix(rnorm(20), ncol=5), name="y")
reps <- rmna(y, ccov=tr, tvcov=dd)
#
# normal intensity, independence model
kalseries(y, intensity="normal", dep="independence", preg=1, pshape=5)
if (FALSE) {
# random effect
kalseries(y, intensity="normal", dep="frailty", preg=1, pinitial=1, psh=5)
# serial dependence
kalseries(y, intensity="normal", dep="serial", preg=1, pinitial=1,
pdep=0.1, psh=5)
# random effect and autoregression
kalseries(y, intensity="normal", dep="frailty", preg=1, pinitial=1,
pdep=0.1, psh=5)
#
# add time-constant variable
kalseries(y, intensity="normal", dep="serial", pinitial=1,
pdep=0.1, psh=5, preg=c(1,0), ccov=treat)
# or equivalently
kalseries(y, intensity="normal", mu=~treat, dep="serial", pinitial=1,
pdep=0.1, psh=5, preg=c(1,0))
# or
kalseries(y, intensity="normal", mu=~b0+b1*treat, dep="serial",
pinitial=1, pdep=0.1, psh=5, preg=c(1,0), envir=reps)
#
# add time-varying variable
kalseries(y, intensity="normal", dep="serial", pinitial=1, pdep=0.1,
psh=5, preg=c(1,0), ccov=treat, ptvc=0, tvc=dose)
# or equivalently, from the environment
dosev <- as.vector(t(dose))
kalseries(y, intensity="normal",
mu=~b0+b1*rep(treat,rep(5,4))+b2*dosev,
dep="serial", pinitial=1, pdep=0.1, psh=5, ptvc=c(1,0,0))
# or from the reps data object
kalseries(y, intensity="normal", mu=~b0+b1*treat+b2*dose,
dep="serial", pinitial=1, pdep=0.1, psh=5, ptvc=c(1,0,0),
envir=reps)
# try power variance family instead of gamma distribution for mixture
kalseries(y, intensity="normal", mu=~b0+b1*treat+b2*dose,
dep="serial", pinitial=1, pdep=0.1, psh=5, ptvc=c(1,0,0),
pfamily=0.1, envir=reps)
# 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.01,1))),
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))
reps <- rmna(conc, ccov=dd, tvcov=tt)
#
# constant shape parameter
kalseries(reps, intensity="gamma", dep="independence", mu=mu,
ptvc=c(-1,-1.1,-1), pshape=1.5)
# or
kalseries(reps, intensity="gamma", dep="independence",
mu=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
ptvc=list(absorption=-1,elimination=-1.1,volume=-1),
pshape=1.2)
# add serial dependence
kalseries(reps, intensity="gamma", dep="serial", pdep=0.9,
mu=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
ptvc=list(absorption=-1,elimination=-1.1,volume=-1),
pshape=0.2)
# time dependent shape parameter
kalseries(reps, intensity="gamma", dep="independence", mu=mu,
shape=shape, ptvc=c(-1,-1.1,-1), pshape=c(-3,0))
# or
kalseries(reps, intensity="gamma", dep="independence",
mu=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
ptvc=list(absorption=-1,elimination=-1.1,volume=-1),
shape=~exp(b1-b2)*times*dose*exp(-exp(b1)*times),
pshape=list(b1=-3,b2=0))
# add serial dependence
kalseries(reps, intensity="gamma", dep="serial", pdep=0.5,
mu=~exp(absorption-volume)*
dose/(exp(absorption)-exp(elimination))*
(exp(-exp(elimination)*times)-exp(-exp(absorption)*times)),
ptvc=list(absorption=-1,elimination=-1.1,volume=-1),
shape=~exp(b1-b2)*times*dose*exp(-exp(b1)*times),
pshape=list(b1=-3,b2=0))
}
Run the code above in your browser using DataLab