| Title: | Non-Normal Repeated Measurements Models |
|---|---|
| Description: | Various functions to fit models for non-normal repeated measurements, such as Binary Random Effects Models with Two Levels of Nesting, Bivariate Beta-binomial Regression Models, Marginal Bivariate Binomial Regression Models, Cormack capture-recapture models, Continuous-time Hidden Markov Chain Models, Discrete-time Hidden Markov Chain Models, Changepoint Location Models using a Continuous-time Two-state Hidden Markov Chain, generalized nonlinear autoregression models, multivariate Gaussian copula models, generalized non-linear mixed models with one random effect, generalized non-linear mixed models using h-likelihood for one random effect, Repeated Measurements Models for Counts with Frailty or Serial Dependence, Repeated Measurements Models for Continuous Variables with Frailty or Serial Dependence, Ordinal Random Effects Models with Dropouts, marginal homogeneity models for square contingency tables, correlated negative binomial models with Kalman update. References include Lindsey's text books, JK Lindsey (2001) <isbn:10-0198508123> and JK Lindsey (1999) <isbn:10-0198505590>. |
| Authors: | Bruce Swihart [cre, aut], Jim Lindsey [aut] (Jim created this package, Bruce is maintaining the CRAN version), T.R. Ten Have [ctb, cph] (Wrote Logit_bin_nest.f90 (in binnest.f).), Richard Cook [ctb, cph] (Wrote calcs.c, calcs.h, calcs1.c, calcs1.h, calcs2.c, calcs2.h, calcs3.c, calcs3.h, calcs4.c, calcs4.h; defs.h; gaps.c, gaps.h.), Iain MacDonald [ctb, cph] (Wrote chidden.f, cphidden.f, hidden.f.), Walter Zucchini [ctb, cph] (Wrote chidden.f, cphidden.f, hidden.f.), Burton Garbow [ctb, cph] (Wrote eigen.f.), Euginia Zharichenko [ctb, cph] (Wrote logitord.f.) |
| Maintainer: | Bruce Swihart <[email protected]> |
| License: | GPL (>=2) |
| Version: | 1.1.10 |
| Built: | 2025-01-01 05:21:10 UTC |
| Source: | https://github.com/swihart/repeated |
binnest is designed to handle binary and binomial data with two
levels of nesting. The first level is the individual and the second will
consist of clusters within individuals.
binnest( response, totals = NULL, nest = NULL, ccov = NULL, tvcov = NULL, mu = ~1, re1 = ~1, re2 = ~1, preg = NULL, pre1 = NULL, pre2 = NULL, binom.mix = c(10, 10), binom.prob = c(0.5, 0.5), fcalls = 900, eps = 0.01, print.level = 0 )binnest( response, totals = NULL, nest = NULL, ccov = NULL, tvcov = NULL, mu = ~1, re1 = ~1, re2 = ~1, preg = NULL, pre1 = NULL, pre2 = NULL, binom.mix = c(10, 10), binom.prob = c(0.5, 0.5), fcalls = 900, eps = 0.01, print.level = 0 )
response |
A list of three column matrices with counts, corresponding
totals (not necessary if the response is binary), and (second-level)
nesting indicator for each individual, one matrix or dataframe of such
counts, or an object of class, response (created by
|
totals |
If |
nest |
If |
ccov |
If |
tvcov |
If |
mu |
If |
re1 |
If |
re2 |
If |
preg |
Initial parameter estimates for the fixed effect regression
model: either the model specified by |
pre1 |
Initial parameter estimates for the first level of nesting
variance model: either the model specified by |
pre2 |
Initial parameter estimates for the second level of nesting
variance model: either the model specified by |
binom.mix |
A vector of two values giving the totals for the binomial distributions used as the mixing distributions at the two levels of nesting. |
binom.prob |
A vector of two values giving the probabilities in the binomial distributions used as the mixing distributions at the two levels of nesting. If they are 0.5, the mixing distributions approximate normal mixing distributions; otherwise, they are skewed. |
fcalls |
Number of function calls allowed. |
eps |
Convergence criterion. |
print.level |
If 1, the iterations are printed out. |
The variance components at the two levels can only depend on the covariates
if response has class, repeated.
A list of classes binnest is returned.
T.R. Ten Have and J.K. Lindsey
Ten Have, T.R., Kunselman, A.R., and Tran, L. (1999) Statistics in Medicine 18, 947-960.
gar, read.list,
restovec, rmna,
tcctomat, tvctomat.
#y <- rbind(matrix(rbinom(20,1,0.6), ncol=4), # matrix(rbinom(20,1,0.4), ncol=4)) y <- matrix(c(1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,0,1,0,1,1,0,1,0, 1,1,1,1,1,1,1,1,0,1,1,0),nrow=10,ncol=4,byrow=TRUE) resp <- restovec(y, nest=1:4, times=FALSE) ccov <- tcctomat(c(rep(0,5),rep(1,5)), name="treatment") reps <- rmna(resp, ccov=ccov) # two random effects binnest(reps, mu=~treatment, preg=c(1,1), pre1=2, pre2=2) # first level random effect only binnest(reps, mu=~treatment, preg=c(1,-1), pre1=1) # second level random effect only binnest(reps, mu=~treatment, preg=c(1,-1), pre2=1)#y <- rbind(matrix(rbinom(20,1,0.6), ncol=4), # matrix(rbinom(20,1,0.4), ncol=4)) y <- matrix(c(1,1,0,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,0,1,0,1,1,0,1,0, 1,1,1,1,1,1,1,1,0,1,1,0),nrow=10,ncol=4,byrow=TRUE) resp <- restovec(y, nest=1:4, times=FALSE) ccov <- tcctomat(c(rep(0,5),rep(1,5)), name="treatment") reps <- rmna(resp, ccov=ccov) # two random effects binnest(reps, mu=~treatment, preg=c(1,1), pre1=2, pre2=2) # first level random effect only binnest(reps, mu=~treatment, preg=c(1,-1), pre1=1) # second level random effect only binnest(reps, mu=~treatment, preg=c(1,-1), pre2=1)
biv.betab fits dependent (logit) linear regression models to a
bivariate beta-binomial distribution.
biv.betab( freq, x = NULL, p, depend = TRUE, print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1 )biv.betab( freq, x = NULL, p, depend = TRUE, print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1 )
freq |
A matrix containing four columns corresponding to 00, 01, 10, and 11 responses. |
x |
A matrix of explanatory variables, containing pairs of columns, one for each response, and the same number of rows as freq. |
p |
Initial parameter estimates: intercept, dependence (if depend is TRUE, and one for each pair of columns of x. |
depend |
If FALSE, the independence (logistic) model is fitted. |
print.level |
Arguments for nlm. |
typsize |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
A list of class bivbetab is returned.
J.K. Lindsey
y <- matrix( c( 2, 1, 1,13, 4, 1, 3, 5, 3, 3, 1, 4, 15, 8, 1, 6),ncol=4,byrow=TRUE) first <- c(0,0,1,1) second <- c(0,1,0,1) self <- cbind(first,second) other <- cbind(second,first) biv.betab(y,cbind(self,other),p=c(-1,2,1,1)) # independence biv.betab(y,cbind(self,other),p=c(-1,1,1),dep=FALSE)y <- matrix( c( 2, 1, 1,13, 4, 1, 3, 5, 3, 3, 1, 4, 15, 8, 1, 6),ncol=4,byrow=TRUE) first <- c(0,0,1,1) second <- c(0,1,0,1) self <- cbind(first,second) other <- cbind(second,first) biv.betab(y,cbind(self,other),p=c(-1,2,1,1)) # independence biv.betab(y,cbind(self,other),p=c(-1,1,1),dep=FALSE)
biv.binom fits (logit) linear regression models to a marginal
bivariate binomial distribution. The covariates must be of length K, that
is the number of 2x2 tables.
biv.binom( freq, marg1 = ~1, marg2 = ~1, interaction = ~1, pmarg1 = 1, pmarg2 = 1, pinteraction = 1, print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1 )biv.binom( freq, marg1 = ~1, marg2 = ~1, interaction = ~1, pmarg1 = 1, pmarg2 = 1, pinteraction = 1, print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1 )
freq |
A four-column matrix containing K 2x2 frequency tables. |
marg1 |
The model formula for the first margin. |
marg2 |
The model formula for the second margin. |
interaction |
The model formula for the interaction. |
pmarg1 |
Initial parameter estimates for the first margin regression. |
pmarg2 |
Initial parameter estimates for the second margin regression. |
pinteraction |
Initial parameter estimates for the interaction regression. |
print.level |
Arguments for nlm. |
typsize |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
A list of class bivbinom is returned.
J.K. Lindsey
# 5 2x2 tables Freq <- matrix(rpois(20,10),ncol=4) x <- c(6,8,10,12,14) print(z <- biv.binom(Freq,marg1=~x,marg2=~x,inter=~x,pmarg1=c(-2,0.08), pmarg2=c(-2,0.1),pinter=c(3,0)))# 5 2x2 tables Freq <- matrix(rpois(20,10),ncol=4) x <- c(6,8,10,12,14) print(z <- biv.binom(Freq,marg1=~x,marg2=~x,inter=~x,pmarg1=c(-2,0.08), pmarg2=c(-2,0.1),pinter=c(3,0)))
capture fits the Cormack capture-recapture model to n sample
periods. Set n to the appropriate value and type eval(setup).
capture(z, n)capture(z, n)
z |
A Poisson generalized linear model object. |
n |
The number of repeated observations. |
n <- periods # number of periods
eval(setup)
This produces the following variables -
p[i]: logit capture probabilities,
pbd: constant capture probability,
d[i]: death parameters,
b[i]: birth parameters,
pw: prior weights.
Then set up a Poisson model for log linear models:
z <- glm(y~model, family=poisson, weights=pw)
and call the function, capture.
If there is constant effort, then all estimates are correct. Otherwise,
n[1], p[1], b[1], are correct only if there is no
birth in period 1. n[s], p[s], are correct only if there is
no death in the last period. phi[s-1] is correct only if effort is
constant in (s-1, s). b[s-1] is correct only if n[s]
and phi[s-1] both are.
capture returns a matrix containing the estimates.
J.K. Lindsey
y <- c(0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,14,1,1,0,2,1,2,1,16,0,2,0,11, 2,13,10,0) n <- 5 eval(setup) # closed population print(z0 <- glm(y~p1+p2+p3+p4+p5, family=poisson, weights=pw)) # deaths and emigration only print(z1 <- update(z0, .~.+d1+d2+d3)) # immigration only print(z2 <- update(z1, .~.-d1-d2-d3+b2+b3+b4)) # deaths, emigration, and immigration print(z3 <- update(z2, .~.+d1+d2+d3)) # add trap dependence print(z4 <- update(z3, .~.+i2+i3)) # constant capture probability over the three middle periods print(z5 <- glm(y~p1+pbd+p5+d1+d2+d3+b2+b3+b4, family=poisson, weights=pw)) # print out estimates capture(z5, n)y <- c(0,1,0,0,0,1,0,1,0,0,0,1,0,0,0,14,1,1,0,2,1,2,1,16,0,2,0,11, 2,13,10,0) n <- 5 eval(setup) # closed population print(z0 <- glm(y~p1+p2+p3+p4+p5, family=poisson, weights=pw)) # deaths and emigration only print(z1 <- update(z0, .~.+d1+d2+d3)) # immigration only print(z2 <- update(z1, .~.-d1-d2-d3+b2+b3+b4)) # deaths, emigration, and immigration print(z3 <- update(z2, .~.+d1+d2+d3)) # add trap dependence print(z4 <- update(z3, .~.+i2+i3)) # constant capture probability over the three middle periods print(z5 <- glm(y~p1+pbd+p5+d1+d2+d3+b2+b3+b4, family=poisson, weights=pw)) # print out estimates capture(z5, n)
catmiss calculates the marginal probabilities of repeated responses.
If there are missing values, it gives both the complete data estimates and
the estimates using all data. It is useful, for example, when a log linear
model is fitted; the resulting fitted values can be supplied to
catmiss to obtain the estimates of the marginal probabilities for
the model. (Note however that the standard errors do not take into account
the fitting of the model.)
catmiss(response, frequency, ccov = NULL)catmiss(response, frequency, ccov = NULL)
response |
A matrix with one column for each of the repeated measures and one row for each possible combination of responses, including the missing values, indicated by NAs. |
frequency |
A vector containing the frequencies. Its length must be a
multiple of the number of rows of |
ccov |
An optional matrix containing the explanatory variables
(time-constant covariates) as columns, with one line per block of responses
in |
A matrix with the probabilities and their standard errors is returned.
J.K. Lindsey
y <- rpois(27,15) r1 <- gl(3,1,27) r2 <- gl(3,3,27) r3 <- gl(3,9) # r1, r2, and r3 are factor variables with 3 indicating missing # independence model with three binary repeated measures # with missing values print(z <- glm(y~r1+r2+r3, family=poisson)) # obtain marginal estimates (no observations with 3 missing values) resp <- cbind(as.integer(r1), as.integer(r2), as.integer(r3))[1:26,] resp <- ifelse(resp==3, NA, resp) catmiss(resp, y[1:26])y <- rpois(27,15) r1 <- gl(3,1,27) r2 <- gl(3,3,27) r3 <- gl(3,9) # r1, r2, and r3 are factor variables with 3 indicating missing # independence model with three binary repeated measures # with missing values print(z <- glm(y~r1+r2+r3, family=poisson)) # obtain marginal estimates (no observations with 3 missing values) resp <- cbind(as.integer(r1), as.integer(r2), as.integer(r3))[1:26,] resp <- ifelse(resp==3, NA, resp) catmiss(resp, y[1:26])
gar fits a first- or second-order generalized autoregression,
possibly with Kalman update over time (first-order only).
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 = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )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 = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )
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,
|
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.) |
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,
|
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, |
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, |
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, |
mu |
A user-specified function of |
shape |
An optional user-specified shape regression function; this may
depend on the location (function) through its second argument, in which
case, |
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, |
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. |
common |
If TRUE, |
preg |
The initial parameter estimates for the location regression
function. If |
pshape |
Zero to two estimates for the shape parameters, depending on
the distribution, if |
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 |
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. |
arch |
If |
transform |
Transformation of the response variable: |
link |
Link function for the mean: |
autocorr |
The form of the (second if two) dependence function:
|
order |
First- or second-order stationary autoregression. |
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
steptol |
Arguments for nlm. |
fscale |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
typsize |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
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 finterp.)
Marginal and individual profiles can be plotted using
mprofile and iprofile and
residuals with plot.residuals.
When the dispersion parameter is not constant over time, volatility
extracts the square root of the dispersion parameter for a fitted model.
A list of classes gar and 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.
J.K. Lindsey
Lindsey, J.K. (1997) Applying Generalized Linear Models. Springer, pp.\ 93–101
Lambert, P. (1996) Statistics in Medicine 15, 1695-1708
# 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) ## Not run: # 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)) ## End(Not run)# 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) ## Not run: # 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)) ## End(Not run)
gausscop fits multivariate repeated measurements models based on the
Gaussian copula with a choice of marginal distributions. Dependence among
responses is provided by the correlation matrix containing random effects
and/or autoregression.
gausscop( response = NULL, distribution = "gamma", mu = NULL, shape = NULL, autocorr = "exponential", pmu = NULL, pshape = NULL, par = NULL, pre = NULL, delta = NULL, shfn = FALSE, common = FALSE, envir = parent.frame(), print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, iterlim = 100, fscale = 1, stepmax = 10 * sqrt(theta %*% theta), typsize = abs(c(theta)) )gausscop( response = NULL, distribution = "gamma", mu = NULL, shape = NULL, autocorr = "exponential", pmu = NULL, pshape = NULL, par = NULL, pre = NULL, delta = NULL, shfn = FALSE, common = FALSE, envir = parent.frame(), print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, iterlim = 100, fscale = 1, stepmax = 10 * sqrt(theta %*% theta), typsize = abs(c(theta)) )
response |
A list of two or three column matrices with response
values, times, and possibly nesting categories, for each individual, one
matrix or dataframe of response values, or an object of class,
|
distribution |
The marginal distribution: exponential, gamma, Weibull, Pareto, inverse Gauss, logistic, Cauchy, Laplace, or Levy. |
mu |
The linear or nonlinear regression model to be fitted for the location parameter. For marginal distributions requiring positive response values, a log link is used. This model can be a function of the parameters or a formula beginning with ~, specifying either a linear regression function for the location parameter in the Wilkinson and Rogers notation or a general function with named unknown parameters that describes the location, returning a vector the same length as the number of observations. |
shape |
The linear or nonlinear regression model to be fitted for the
log shape parameter. This can be a function of the parameters or a formula
beginning with ~, specifying either a linear regression function for the
location parameter in the Wilkinson and Rogers notation or a general
function with named unknown parameters that describes the location. If it
contains unknown parameters, the keyword |
autocorr |
The form of the autocorrelation function:
|
pmu |
Initial parameter estimates for the location regression model. |
pshape |
Initial parameter estimate for the shape regression model. |
par |
If supplied, an initial estimate for the autocorrelation parameter. |
pre |
Zero, one or two parameter estimates for the variance components, depending on the number of levels of nesting. |
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, |
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. |
common |
If TRUE, |
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
typsize |
Arguments for nlm. |
With two levels of nesting, the first is the individual and the second will consist of clusters within individuals.
For clustered (non-longitudinal) data, where only random effects will be
fitted, times are not necessary.
This function 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 in which case factor variables cannot be used and parameters
must be scalars. (See finterp.)
A list of class gausscop is returned that contains all of
the relevant information calculated, including error codes.
J.K. Lindsey
Song, P.X.K. (2000) Multivariate dispersion models generated from Gaussian copula. Scandinavian Journal of Statistics 27, 305-320.
# linear models y <- matrix(rgamma(40,1,1),ncol=5)+rep(rgamma(8,0.5,1),5) x1 <- c(rep(0,4),rep(1,4)) reps <- rmna(restovec(y),ccov=tcctomat(x1)) # independence with default gamma marginals # compare with gnlm::gnlr(y, pmu=1, psh=0, dist="gamma", env=reps) gausscop(y, pmu=1, pshape=0, env=reps) gausscop(y, mu=~x1, pmu=c(1,0), pshape=0, env=reps) # AR(1) gausscop(y, pmu=1, pshape=0, par=0.1, env=reps) ## Not run: # random effect gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps) # try other marginal distributions gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps, dist="Weibull") gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps, dist="inverse Gauss", stepmax=1) gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps, dist="Cauchy") # # first-order one-compartment model # create 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) lmu <- function(p) p[1]-p[3]+log(dose/(exp(p[1])-exp(p[2]))* (exp(-exp(p[2])*times)-exp(-exp(p[1])*times))) lshape <- function(p) p[1]-p[2]+log(times*dose)-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") conc <- matrix(c(3.65586845,0.01000000,0.01000000,0.01731192,1.68707608, 0.01000000,4.67338974,4.79679942,1.86429851,1.82886732,1.54708795, 0.57592054,0.08014232,0.09436425,0.26106139,0.11125534,0.22685364, 0.22896015,0.04886441,0.01000000,33.59011263,16.89115866,19.99638316, 16.94021361,9.95440037,7.10473948,2.97769676,1.53785279,2.13059515, 0.72562344,1.27832563,1.33917155,0.99811111,0.23437424,0.42751355, 0.65702300,0.41126684,0.15406463,0.03092312,0.14672610), ncol=20,byrow=TRUE) conc <- restovec(conc) reps <- rmna(conc, ccov=dd, tvcov=tt) # constant shape parameter gausscop(conc, mu=lmu, pmu=log(c(1,0.4,0.1)), par=0.5, pshape=0, envir=reps) # or gausscop(conc, mu=~absorption-volume+ log(dose/(exp(absorption)-exp(elimination))* (exp(-exp(elimination)*times)-exp(-exp(absorption)*times))), pmu=list(absorption=0,elimination=log(0.4),volume=log(0.1)), par=0.5, pshape=0, envir=reps) # compare to gar autoregression gar(conc, dist="gamma", times=1:20, mu=mu, preg=log(c(1,0.4,0.1)), pdepend=0.5, pshape=1) # # time dependent shape parameter gausscop(conc, mu=lmu, shape=lshape, pmu=log(c(1,0.4,0.1)), par=0.5, pshape=c(-0.1,-0.1)) # or gausscop(conc, mu=~absorption-volume+ log(dose/(exp(absorption)-exp(elimination))* (exp(-exp(elimination)*times)-exp(-exp(absorption)*times))), shape=~b1-b2+log(times*dose)-exp(b1)*times, pmu=list(absorption=0,elimination=log(0.4),volume=log(0.1)), par=0.5, pshape=list(b1=-0.1,b2=-0.1), envir=reps) # # shape depends on location lshape <- function(p, mu) p[1]*log(abs(mu)) gausscop(conc, mu=lmu, shape=lshape, shfn=TRUE, pmu=log(c(1,0.4,0.1)), par=0.5, pshape=1) # or gausscop(conc, mu=~absorption-volume+ log(dose/(exp(absorption)-exp(elimination))* (exp(-exp(elimination)*times)-exp(-exp(absorption)*times))), shape=~d*log(abs(mu)), shfn=TRUE, pmu=list(absorption=0,elimination=log(0.4),volume=log(0.1)), par=0.5, pshape=list(d=1), envir=reps) ## End(Not run)# linear models y <- matrix(rgamma(40,1,1),ncol=5)+rep(rgamma(8,0.5,1),5) x1 <- c(rep(0,4),rep(1,4)) reps <- rmna(restovec(y),ccov=tcctomat(x1)) # independence with default gamma marginals # compare with gnlm::gnlr(y, pmu=1, psh=0, dist="gamma", env=reps) gausscop(y, pmu=1, pshape=0, env=reps) gausscop(y, mu=~x1, pmu=c(1,0), pshape=0, env=reps) # AR(1) gausscop(y, pmu=1, pshape=0, par=0.1, env=reps) ## Not run: # random effect gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps) # try other marginal distributions gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps, dist="Weibull") gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps, dist="inverse Gauss", stepmax=1) gausscop(y, pmu=1, pshape=0, pre=0.1, env=reps, dist="Cauchy") # # first-order one-compartment model # create 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) lmu <- function(p) p[1]-p[3]+log(dose/(exp(p[1])-exp(p[2]))* (exp(-exp(p[2])*times)-exp(-exp(p[1])*times))) lshape <- function(p) p[1]-p[2]+log(times*dose)-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") conc <- matrix(c(3.65586845,0.01000000,0.01000000,0.01731192,1.68707608, 0.01000000,4.67338974,4.79679942,1.86429851,1.82886732,1.54708795, 0.57592054,0.08014232,0.09436425,0.26106139,0.11125534,0.22685364, 0.22896015,0.04886441,0.01000000,33.59011263,16.89115866,19.99638316, 16.94021361,9.95440037,7.10473948,2.97769676,1.53785279,2.13059515, 0.72562344,1.27832563,1.33917155,0.99811111,0.23437424,0.42751355, 0.65702300,0.41126684,0.15406463,0.03092312,0.14672610), ncol=20,byrow=TRUE) conc <- restovec(conc) reps <- rmna(conc, ccov=dd, tvcov=tt) # constant shape parameter gausscop(conc, mu=lmu, pmu=log(c(1,0.4,0.1)), par=0.5, pshape=0, envir=reps) # or gausscop(conc, mu=~absorption-volume+ log(dose/(exp(absorption)-exp(elimination))* (exp(-exp(elimination)*times)-exp(-exp(absorption)*times))), pmu=list(absorption=0,elimination=log(0.4),volume=log(0.1)), par=0.5, pshape=0, envir=reps) # compare to gar autoregression gar(conc, dist="gamma", times=1:20, mu=mu, preg=log(c(1,0.4,0.1)), pdepend=0.5, pshape=1) # # time dependent shape parameter gausscop(conc, mu=lmu, shape=lshape, pmu=log(c(1,0.4,0.1)), par=0.5, pshape=c(-0.1,-0.1)) # or gausscop(conc, mu=~absorption-volume+ log(dose/(exp(absorption)-exp(elimination))* (exp(-exp(elimination)*times)-exp(-exp(absorption)*times))), shape=~b1-b2+log(times*dose)-exp(b1)*times, pmu=list(absorption=0,elimination=log(0.4),volume=log(0.1)), par=0.5, pshape=list(b1=-0.1,b2=-0.1), envir=reps) # # shape depends on location lshape <- function(p, mu) p[1]*log(abs(mu)) gausscop(conc, mu=lmu, shape=lshape, shfn=TRUE, pmu=log(c(1,0.4,0.1)), par=0.5, pshape=1) # or gausscop(conc, mu=~absorption-volume+ log(dose/(exp(absorption)-exp(elimination))* (exp(-exp(elimination)*times)-exp(-exp(absorption)*times))), shape=~d*log(abs(mu)), shfn=TRUE, pmu=list(absorption=0,elimination=log(0.4),volume=log(0.1)), par=0.5, pshape=list(d=1), envir=reps) ## End(Not run)
glmm fits a generalized linear mixed model with a random intercept
using a normal mixing distribution computed by Gauss-Hermite integration.
For the normal, gamma, and inverse Gaussian distributions, the deviances
supplied are -2 log likelihood, not the usual glm deviance;
the degrees of freedom take into account estimation of the dispersion
parameter.
glmm( formula, family = gaussian, data = list(), weights = NULL, offset = NULL, nest, delta = 1, maxiter = 20, points = 10, print.level = 0, control = glm.control(epsilon = 1e-04, maxit = 10, trace = FALSE) )glmm( formula, family = gaussian, data = list(), weights = NULL, offset = NULL, nest, delta = 1, maxiter = 20, points = 10, print.level = 0, control = glm.control(epsilon = 1e-04, maxit = 10, trace = FALSE) )
formula |
A symbolic description of the model to be fitted. If it contains transformations of the data, including cbind for binomial data, a dataframe must be supplied. |
family |
A description of the error distribution and link function to
be used in the model; see |
data |
A dataframe containing the variables in the model, that is optional in simple cases, but required in certain situations as specified elsewhere in this help page. |
weights |
An optional weight vector. If this is used, data must be supplied in a data.frame. |
offset |
The known component in the linear predictor. If this is used, data must be supplied in a data.frame. An offset cannot be specified in the model formula. |
nest |
The variable classifying observations by the unit (cluster) upon which they were observed. |
delta |
If the response variable has been transformed, this is the Jacobian of that transformation, so that AICs are comparable. |
maxiter |
The maximum number of iterations of the outer loop for numerical integration. |
points |
The number of points for Gauss-Hermite integration of the random effect. |
print.level |
If set equal to 2, the log probabilities are printed out when the underflow error is given. |
control |
A list of parameters for controlling the fitting process. |
If weights and/or offset are to be used or the formula transforms some
variables, all of the data must be supplied in a dataframe. Because the
glm function is such a hack, if this is not done, weird error
messages will result.
na.omit is not allowed.
glmm returns a list of class glmm
J.K. Lindsey
# Poisson counts nest <- gl(5,4) y <- rpois(20,5+2*as.integer(nest)) # overdispersion model glmm(y~1, family=poisson, nest=gl(20,1), points=3) # clustered model glmm(y~1, family=poisson, nest=nest, points=3) # # binomial data with model for overdispersion df <- data.frame(r=rbinom(10,10,0.5), n=rep(10,10), x=c(rep(0,5), rep(1,5)), nest=1:10) glmm(cbind(r,n-r)~x, family=binomial, nest=nest, data=df)# Poisson counts nest <- gl(5,4) y <- rpois(20,5+2*as.integer(nest)) # overdispersion model glmm(y~1, family=poisson, nest=gl(20,1), points=3) # clustered model glmm(y~1, family=poisson, nest=nest, points=3) # # binomial data with model for overdispersion df <- data.frame(r=rbinom(10,10,0.5), n=rep(10,10), x=c(rep(0,5), rep(1,5)), nest=1:10) glmm(cbind(r,n-r)~x, family=binomial, nest=nest, data=df)
gnlmix fits user-specified nonlinear regression equations to one or
both parameters of the common one and two parameter distributions. One
parameter of the location regression is random with some specified mixing
distribution.
gnlmix( y = NULL, distribution = "normal", mixture = "normal", random = NULL, nest = NULL, mu = NULL, shape = NULL, linear = NULL, pmu = NULL, pshape = NULL, pmix = NULL, delta = 1, common = FALSE, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1, eps = 1e-04, points = 5, steps = 10 )gnlmix( y = NULL, distribution = "normal", mixture = "normal", random = NULL, nest = NULL, mu = NULL, shape = NULL, linear = NULL, pmu = NULL, pshape = NULL, pmix = NULL, delta = 1, common = FALSE, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1, eps = 1e-04, points = 5, steps = 10 )
y |
A response vector of uncensored data, a two column matrix for
binomial data, or an object of class, |
distribution |
The distribution for the response: binomial, beta binomial, double binomial, mult(iplicative) binomial, Poisson, negative binomial, double Poisson, mult(iplicative) Poisson, gamma count, Consul generalized Poisson, logarithmic series, geometric, normal, inverse Gauss, logistic, exponential, gamma, Weibull, extreme value, Cauchy, Pareto, Laplace, Levy, beta, simplex, or two-sided power. (For definitions of distributions, see the corresponding [dpqr]distribution help.) |
mixture |
The mixing distribution for the random parameter: normal, Cauchy, logistic, Laplace, inverse Gauss, gamma, inverse gamma, Weibull, beta, simplex, or two-sided power. The first four have zero location parameter, the next three have unit location parameter, and the last two have location parameter set to 0.5. |
random |
The name of the random parameter in the |
nest |
The variable classifying observations by the unit upon which
they were observed. Ignored if |
mu |
A user-specified formula containing named unknown parameters,
giving the regression equation for the location parameter. This may contain
the keyword, |
shape |
A user-specified formula containing named unknown parameters,
giving the regression equation for the shape parameter. This may contain
the keyword, |
linear |
A formula beginning with ~ in W&R notation, specifying the linear part of the regression function for the location parameter or list of two such expressions for the location and/or shape parameters. |
pmu |
Vector of initial estimates for the location parameters. These must be supplied either in their order of appearance in the formula or in a named list. |
pshape |
Vector of initial estimates for the shape parameters. These must be supplied either in their order of appearance in the expression or in a named list. |
pmix |
Initial estimate for the logarithm of the dispersion parameter of the mixing distribution. |
delta |
Scalar or vector giving the unit of measurement (always one
for discrete data) for each response value, set to unity by default. For
example, if a response is measured to two decimals, |
common |
If TRUE, the formulae with unknowns for the location and
shape have names in common. All parameter estimates must be supplied in
|
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
typsize |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
eps |
Precision of the Romberg integration. |
points |
For the Romberg integration, the number of extrapolation points so that 2*points is the order of integration, by default set to 5; points=2 is Simpson's rule. |
steps |
For the Romberg integration, the maximum number of steps, by default set to 10. |
It is recommended that initial estimates for pmu and pshape
be obtained from gnlr.
These nonlinear regression models must be supplied as formulae where
parameters are unknowns. (See finterp.)
A list of class gnlm is returned that contains all of the
relevant information calculated, including error codes.
J.K. Lindsey
dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8) #y <- rgamma(20,shape=2+0.3*dose,scale=2)+rep(rnorm(4,0,4),rep(5,4)) y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699, 4.359469, 1.900681, 17.425948, 4.503345, 2.691792, 5.731100, 10.534971, 11.220260, 6.968932, 4.094357, 16.393806, 14.656584, 8.786133, 20.972267, 17.178012) resp <- restovec(matrix(y, nrow=4, byrow=TRUE), name="y") reps <- rmna(resp, tvcov=tvctomat(matrix(dose, nrow=4, byrow=TRUE), name="dose")) # same linear normal model with random normal intercept fitted four ways # compare with growth::elliptic(reps, model=~dose, preg=c(0,0.6), pre=4) glmm(y~dose, nest=individuals, data=reps) gnlmm(reps, mu=~dose, pmu=c(8.7,0.25), psh=3.5, psd=3) gnlmix(reps, mu=~a+b*dose+rand, random="rand", pmu=c(8.7,0.25), pshape=3.44, pmix=2.3) ## Not run: # gamma model with log link and random normal intercept fitted three ways glmm(y~dose, family=Gamma(link=log), nest=individuals, data=reps, points=8) gnlmm(reps, distribution="gamma", mu=~exp(a+b*dose), pmu=c(2,0.03), psh=1, psd=0.3) gnlmix(reps, distribution="gamma", mu=~exp(a+b*dose+rand), random="rand", pmu=c(2,0.04), pshape=1, pmix=-2) # gamma model with log link and random gamma mixtures gnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a*rand+b*dose), random="rand", pmu=c(2,0.04), pshape=1.24, pmix=3.5) gnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a+b*dose)*rand, random="rand", pmu=c(2,0.04), pshape=1.24, pmix=2.5) ## End(Not run)dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8) #y <- rgamma(20,shape=2+0.3*dose,scale=2)+rep(rnorm(4,0,4),rep(5,4)) y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699, 4.359469, 1.900681, 17.425948, 4.503345, 2.691792, 5.731100, 10.534971, 11.220260, 6.968932, 4.094357, 16.393806, 14.656584, 8.786133, 20.972267, 17.178012) resp <- restovec(matrix(y, nrow=4, byrow=TRUE), name="y") reps <- rmna(resp, tvcov=tvctomat(matrix(dose, nrow=4, byrow=TRUE), name="dose")) # same linear normal model with random normal intercept fitted four ways # compare with growth::elliptic(reps, model=~dose, preg=c(0,0.6), pre=4) glmm(y~dose, nest=individuals, data=reps) gnlmm(reps, mu=~dose, pmu=c(8.7,0.25), psh=3.5, psd=3) gnlmix(reps, mu=~a+b*dose+rand, random="rand", pmu=c(8.7,0.25), pshape=3.44, pmix=2.3) ## Not run: # gamma model with log link and random normal intercept fitted three ways glmm(y~dose, family=Gamma(link=log), nest=individuals, data=reps, points=8) gnlmm(reps, distribution="gamma", mu=~exp(a+b*dose), pmu=c(2,0.03), psh=1, psd=0.3) gnlmix(reps, distribution="gamma", mu=~exp(a+b*dose+rand), random="rand", pmu=c(2,0.04), pshape=1, pmix=-2) # gamma model with log link and random gamma mixtures gnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a*rand+b*dose), random="rand", pmu=c(2,0.04), pshape=1.24, pmix=3.5) gnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a+b*dose)*rand, random="rand", pmu=c(2,0.04), pshape=1.24, pmix=2.5) ## End(Not run)
gnlmm fits user-specified nonlinear regression equations to one or
both parameters of the common one and two parameter distributions. The
intercept of the location regression has a normally-distributed random
effect. This normal mixing distribution is computed by Gauss-Hermite
integration.
gnlmm( y = NULL, distribution = "normal", mu = NULL, shape = NULL, linear = NULL, nest = NULL, pmu = NULL, pshape = NULL, psd = NULL, exact = FALSE, wt = 1, delta = 1, shfn = FALSE, scale = NULL, points = 10, common = FALSE, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = sqrt(p %*% p)/10, steptol = 1e-05, iterlim = 100, fscale = 1 )gnlmm( y = NULL, distribution = "normal", mu = NULL, shape = NULL, linear = NULL, nest = NULL, pmu = NULL, pshape = NULL, psd = NULL, exact = FALSE, wt = 1, delta = 1, shfn = FALSE, scale = NULL, points = 10, common = FALSE, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = sqrt(p %*% p)/10, steptol = 1e-05, iterlim = 100, fscale = 1 )
y |
A response vector for uncensored data, a two column matrix for
binomial data or censored data, with the second column being the censoring
indicator (1: uncensored, 0: right censored, -1: left censored), or an
object of class, |
distribution |
Either a character string containing the name of the distribution or a function giving the -log likelihood and calling the location and shape functions. Distributions are binomial, beta binomial, double binomial, mult(iplicative) binomial, Poisson, negative binomial, double Poisson, mult(iplicative) Poisson, gamma count, Consul generalized Poisson, logarithmic series, geometric, normal, inverse Gauss, logistic, exponential, gamma, Weibull, extreme value, Cauchy, Pareto, Laplace, and Levy, beta, simplex, and two-sided power. All but the binomial-based distributions and the beta, simplex, and two-sided power may be right and/or left censored. (For definitions of distributions, see the corresponding [dpqr]distribution help.) |
mu |
A user-specified function of |
shape |
A user-specified function of |
linear |
A formula beginning with ~ in W&R notation, specifying the linear part of the regression function for the location parameter or list of two such expressions for the location and/or shape parameters. |
nest |
The variable classifying observations by the unit upon which
they were observed. Ignored if |
pmu |
Vector of initial estimates for the location parameters. If
|
pshape |
Vector of initial estimates for the shape parameters. If
|
psd |
Initial estimate of the standard deviation of the normal mixing distribution. |
exact |
If TRUE, fits the exact likelihood function for continuous data by integration over intervals of observation, i.e. interval censoring. |
wt |
Weight vector. |
delta |
Scalar or vector giving the unit of measurement (always one
for discrete data) for each response value, set to unity by default.
Ignored if y has class, response. For example, if a response is measured to
two decimals, |
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. |
scale |
The scale on which the random effect is applied:
|
points |
The number of points for Gauss-Hermite integration of the random effect. |
common |
If TRUE, |
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
typsize |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
The scale of the random effect is the link function to be applied.
For example, if it is log, the supplied mean function, mu, is
transformed as exp(log(mu)+sd), where sd is the random effect parameter.
It is recommended that initial estimates for pmu and pshape
be obtained from gnlr.
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 finterp.)
The printed output includes the -log likelihood (not the deviance), the corresponding AIC, the maximum likelihood estimates, standard errors, and correlations.
A list of class gnlm is returned that contains all of the
relevant information calculated, including error codes.
J.K. Lindsey
# data objects sex <- c(0,1,1) sx <- tcctomat(sex) dose <- matrix(rpois(30,10),nrow=3) dd <- tvctomat(dose) # vectors for functions dose <- as.vector(t(dose)) sex <- c(rep(0,10),rep(1,20)) nest <- rbind(rep(1,10),rep(2,10),rep(3,10)) #y <- rgamma(30,2,scale=exp(0.2+0.1*dose+0.1*sex+rep(rnorm(3),rep(10,3)))/2) y <- c(0.6490851,0.9313931,0.4765569,0.4188045,2.8339637,2.8158090, 2.6059975,2.9958184,2.7351583,3.2884980,1.1180961,0.9443986,1.7915571, 9.0013379,2.3969570,3.4227356,0.5045518,0.7452521,1.8712467,3.6814198, 0.1489849,1.0327552,0.6102406,1.1536620,2.9145237,9.2847798,5.6454605, 1.9759672,1.5798008,5.1024496) y <- restovec(matrix(y, nrow=3), nest=nest, name="y") reps <- rmna(y, ccov=sx, tvcov=dd) # # log linear regression with gamma distribution mu <- function(p) exp(p[1]+p[2]*sex+p[3]*dose) ## print(z <- gnlm::gnlr(y, dist="gamma", mu=mu, pmu=c(1,0,0), pshape=1)) ## starting values for pmu and pshape from z$coef[1:3] and z$coef[4] respectively gnlmm(y, dist="gamma", mu=mu, nest=nest, pmu=10*c(0.59072535, 0.32618702, 0.01024245),pshape=1, psd=0.1, points=3) # or equivalently gnlmm(y, dist="gamma", mu=~exp(b0+b1*sex+b2*dose), nest=nest, pmu=10*c(0.59072535, 0.32618702, 0.01024245),pshape=1, psd=0.1, points=3, envir=reps) ## Not run: # or with identity link print(z <- gnlm::gnlr(y, dist="gamma", mu=~sex+dose, pmu=c(0.1,0,0), pshape=1)) gnlmm(y, dist="gamma", mu=~sex+dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=~b0+b1*sex+b2*dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], psd=0.1, points=3, envir=reps) # # nonlinear regression with gamma distribution mu <- function(p) p[1]+exp(p[2]+p[3]*sex+p[4]*dose) print(z <- gnlm::gnlr(y, dist="gamma", mu=mu, pmu=c(1,1,0,0), pshape=1)) gnlmm(y, dist="gamma", mu=mu, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], psd=0.1, points=3) # or mu2 <- function(p, linear) p[1]+exp(linear) gnlmm(y, dist="gamma", mu=mu2, linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=1, psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=~a+exp(linear), linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=1, psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=~b4+exp(b0+b1*sex+b2*dose), nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], psd=0.1, points=3, envir=reps) # # include regression for the shape parameter with same mu function shape <- function(p) p[1]+p[2]*sex print(z <- gnlm::gnlr(y, dist="gamma", mu=mu, shape=shape, pmu=z$coef[1:4], pshape=rep(1,2))) gnlmm(y, dist="gamma", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], psd=0.1, points=3, envir=reps) # or gnlmm(y, dist="gamma", mu=~b4+exp(b0+b1*sex+b2*dose), shape=~a1+a2*sex, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], psd=0.1, points=3, envir=reps) ## End(Not run)# data objects sex <- c(0,1,1) sx <- tcctomat(sex) dose <- matrix(rpois(30,10),nrow=3) dd <- tvctomat(dose) # vectors for functions dose <- as.vector(t(dose)) sex <- c(rep(0,10),rep(1,20)) nest <- rbind(rep(1,10),rep(2,10),rep(3,10)) #y <- rgamma(30,2,scale=exp(0.2+0.1*dose+0.1*sex+rep(rnorm(3),rep(10,3)))/2) y <- c(0.6490851,0.9313931,0.4765569,0.4188045,2.8339637,2.8158090, 2.6059975,2.9958184,2.7351583,3.2884980,1.1180961,0.9443986,1.7915571, 9.0013379,2.3969570,3.4227356,0.5045518,0.7452521,1.8712467,3.6814198, 0.1489849,1.0327552,0.6102406,1.1536620,2.9145237,9.2847798,5.6454605, 1.9759672,1.5798008,5.1024496) y <- restovec(matrix(y, nrow=3), nest=nest, name="y") reps <- rmna(y, ccov=sx, tvcov=dd) # # log linear regression with gamma distribution mu <- function(p) exp(p[1]+p[2]*sex+p[3]*dose) ## print(z <- gnlm::gnlr(y, dist="gamma", mu=mu, pmu=c(1,0,0), pshape=1)) ## starting values for pmu and pshape from z$coef[1:3] and z$coef[4] respectively gnlmm(y, dist="gamma", mu=mu, nest=nest, pmu=10*c(0.59072535, 0.32618702, 0.01024245),pshape=1, psd=0.1, points=3) # or equivalently gnlmm(y, dist="gamma", mu=~exp(b0+b1*sex+b2*dose), nest=nest, pmu=10*c(0.59072535, 0.32618702, 0.01024245),pshape=1, psd=0.1, points=3, envir=reps) ## Not run: # or with identity link print(z <- gnlm::gnlr(y, dist="gamma", mu=~sex+dose, pmu=c(0.1,0,0), pshape=1)) gnlmm(y, dist="gamma", mu=~sex+dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=~b0+b1*sex+b2*dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], psd=0.1, points=3, envir=reps) # # nonlinear regression with gamma distribution mu <- function(p) p[1]+exp(p[2]+p[3]*sex+p[4]*dose) print(z <- gnlm::gnlr(y, dist="gamma", mu=mu, pmu=c(1,1,0,0), pshape=1)) gnlmm(y, dist="gamma", mu=mu, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], psd=0.1, points=3) # or mu2 <- function(p, linear) p[1]+exp(linear) gnlmm(y, dist="gamma", mu=mu2, linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=1, psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=~a+exp(linear), linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=1, psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=~b4+exp(b0+b1*sex+b2*dose), nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], psd=0.1, points=3, envir=reps) # # include regression for the shape parameter with same mu function shape <- function(p) p[1]+p[2]*sex print(z <- gnlm::gnlr(y, dist="gamma", mu=mu, shape=shape, pmu=z$coef[1:4], pshape=rep(1,2))) gnlmm(y, dist="gamma", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], psd=0.1, points=3) # or gnlmm(y, dist="gamma", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], psd=0.1, points=3, envir=reps) # or gnlmm(y, dist="gamma", mu=~b4+exp(b0+b1*sex+b2*dose), shape=~a1+a2*sex, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], psd=0.1, points=3, envir=reps) ## End(Not run)
gnlmm3 fits user-specified nonlinear regression equations to one or
more parameters of the common three parameter distributions. The intercept
of the location regression has a normally-distributed random effect. This
normal mixing distribution is computed by Gauss-Hermite integration.
gnlmm3( y = NULL, distribution = "normal", mu = NULL, shape = NULL, nest = NULL, family = NULL, linear = NULL, pmu = NULL, pshape = NULL, pfamily = NULL, psd = NULL, exact = FALSE, wt = 1, scale = NULL, points = 10, common = FALSE, delta = 1, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1 )gnlmm3( y = NULL, distribution = "normal", mu = NULL, shape = NULL, nest = NULL, family = NULL, linear = NULL, pmu = NULL, pshape = NULL, pfamily = NULL, psd = NULL, exact = FALSE, wt = 1, scale = NULL, points = 10, common = FALSE, delta = 1, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1 )
y |
A response vector for uncensored data, a two column matrix for
binomial data or censored data, with the second column being the censoring
indicator (1: uncensored, 0: right censored, -1: left censored), or an
object of class, |
distribution |
Either a character string containing the name of the distribution or a function giving the -log likelihood and calling the location, shape, and family functions. Distributions are Box-Cox transformed normal, generalized inverse Gauss, generalized logistic, Hjorth, generalized gamma, Burr, generalized Weibull, power exponential, Student t, generalized extreme value, power variance function Poisson, and skew Laplace. (For definitions of distributions, see the corresponding [dpqr]distribution help.) |
mu |
A user-specified function of |
shape |
A user-specified function of |
nest |
The variable classifying observations by the unit upon which
they were observed. Ignored if |
family |
A user-specified function of |
linear |
A formula beginning with ~ in W&R notation, specifying the linear part of the regression function for the location parameter or list of two such expressions for the location and/or shape parameters. |
pmu |
Vector of initial estimates for the location parameters. If
|
pshape |
Vector of initial estimates for the shape parameters. If
|
pfamily |
Vector of initial estimates for the family parameters. If
|
psd |
Initial estimate of the standard deviation of the normal mixing distribution. |
exact |
If TRUE, fits the exact likelihood function for continuous data by integration over intervals of observation, i.e. interval censoring. |
wt |
Weight vector. |
scale |
The scale on which the random effect is applied:
|
points |
The number of points for Gauss-Hermite integration of the random effect. |
common |
If TRUE, at least two of |
delta |
Scalar or vector giving the unit of measurement (always one
for discrete data) for each response value, set to unity by default.
Ignored if y has class, response. For example, if a response is measured to
two decimals, |
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
typsize |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
The scale of the random effect is the link function to be applied.
For example, if it is log, the supplied mean function, mu, is
transformed as exp(log(mu)+sd), where sd is the random effect parameter.
It is recommended that initial estimates for pmu, pshape, and
pfamily be obtained from gnlr3.
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 finterp.)
The printed output includes the -log likelihood (not the deviance), the corresponding AIC, the maximum likelihood estimates, standard errors, and correlations.
A list of class gnlm is returned that contains all of the
relevant information calculated, including error codes.
J.K. Lindsey
# data objects sex <- c(0,1,1) sx <- tcctomat(sex) #dose <- matrix(rpois(30,10),nrow=3) dose <- matrix(c(8,9,11,9,11,11,7,8,7,12,8,8,9,10,15,10,9,9,20,14,4,7, 4,13,10,13,6,13,11,17),nrow=3) dd <- tvctomat(dose) # vectors for functions dose <- as.vector(t(dose)) sex <- c(rep(0,10),rep(1,20)) nest <- rbind(rep(1,10),rep(2,10),rep(3,10)) #y <- (rt(30,5)+exp(0.2+0.3*dose+0.5*sex+rep(rnorm(3),rep(10,3))))*3 y <- c(62.39712552,196.94419614,2224.74940087,269.56691601,12.86079662, 14.96743546, 47.45765042,156.51381687,508.68804438,281.11065302, 92.32443655, 81.88000484, 40.26357733, 13.04433670, 15.58490237, 63.62154867, 23.69677549, 53.52885894, 88.02507682, 34.04302506, 44.28232323,116.80732423,106.72564484, 25.09749055, 12.61839145, -0.04060996,153.32670123, 63.25866087, 17.79852591,930.52558064) y <- restovec(matrix(y, nrow=3), nest=nest, name="y") reps <- rmna(y, ccov=sx, tvcov=dd) # # log linear regression with Student t distribution mu <- function(p) exp(p[1]+p[2]*sex+p[3]*dose) ## print(z <- gnlm::gnlr3(y, dist="Student", mu=mu, pmu=c(0,0,0), pshape=1, pfamily=1)) ## starting values for pmu and pshape from z$coef[1:3] and z$coef[4] respectively ## starting value for pfamily in z$coef[5] gnlmm3(y, dist="Student", mu=mu, nest=nest, pmu=c(3.69,-1.19, 0.039), pshape=5, pfamily=0, psd=40, points=3) # or equivalently gnlmm3(y, dist="Student", mu=~exp(b0+b1*sex+b2*dose), nest=nest, pmu=c(3.69,-1.19, 0.039), pshape=5, pfamily=0, psd=40, points=3, envir=reps) ## Not run: # or with identity link print(z <- gnlm::gnlr3(y, dist="Student", mu=~sex+dose, pmu=c(0.1,0,0), pshape=1, pfamily=1)) gnlmm3(y, dist="Student", mu=~sex+dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], pfamily=z$coef[5], psd=50, points=3) # or gnlmm3(y, dist="Student", mu=~b0+b1*sex+b2*dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], pfamily=z$coef[5], psd=50, points=3, envir=reps) # # nonlinear regression with Student t distribution mu <- function(p) p[1]+exp(p[2]+p[3]*sex+p[4]*dose) print(z <- gnlm::gnlr3(y, dist="Student", mu=mu, pmu=c(1,1,0,0), pshape=1, pfamily=1)) gnlmm3(y, dist="Student", mu=mu, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3) # or mu2 <- function(p, linear) p[1]+exp(linear) gnlmm3(y, dist="Student", mu=mu2, linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3) # or gnlmm3(y, dist="Student", mu=~a+exp(linear), linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3) # or gnlmm3(y, dist="Student", mu=~b4+exp(b0+b1*sex+b2*dose), nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3, envir=reps) # # include regression for the shape parameter with same mu function shape <- function(p) p[1]+p[2]*sex print(z <- gnlm::gnlr3(y, dist="Student", mu=mu, shape=shape, pmu=z$coef[1:4], pshape=c(z$coef[5],0), pfamily=z$coef[6])) gnlmm3(y, dist="Student", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], pfamily=z$coef[7], psd=5, points=3) # or gnlmm3(y, dist="Student", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], pfamily=z$coef[7], psd=5, points=3, envir=reps) # or gnlmm3(y, dist="Student", mu=~b4+exp(b0+b1*sex+b2*dose), shape=~a1+a2*sex, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], pfamily=z$coef[7], psd=5, points=3, envir=reps) ## End(Not run)# data objects sex <- c(0,1,1) sx <- tcctomat(sex) #dose <- matrix(rpois(30,10),nrow=3) dose <- matrix(c(8,9,11,9,11,11,7,8,7,12,8,8,9,10,15,10,9,9,20,14,4,7, 4,13,10,13,6,13,11,17),nrow=3) dd <- tvctomat(dose) # vectors for functions dose <- as.vector(t(dose)) sex <- c(rep(0,10),rep(1,20)) nest <- rbind(rep(1,10),rep(2,10),rep(3,10)) #y <- (rt(30,5)+exp(0.2+0.3*dose+0.5*sex+rep(rnorm(3),rep(10,3))))*3 y <- c(62.39712552,196.94419614,2224.74940087,269.56691601,12.86079662, 14.96743546, 47.45765042,156.51381687,508.68804438,281.11065302, 92.32443655, 81.88000484, 40.26357733, 13.04433670, 15.58490237, 63.62154867, 23.69677549, 53.52885894, 88.02507682, 34.04302506, 44.28232323,116.80732423,106.72564484, 25.09749055, 12.61839145, -0.04060996,153.32670123, 63.25866087, 17.79852591,930.52558064) y <- restovec(matrix(y, nrow=3), nest=nest, name="y") reps <- rmna(y, ccov=sx, tvcov=dd) # # log linear regression with Student t distribution mu <- function(p) exp(p[1]+p[2]*sex+p[3]*dose) ## print(z <- gnlm::gnlr3(y, dist="Student", mu=mu, pmu=c(0,0,0), pshape=1, pfamily=1)) ## starting values for pmu and pshape from z$coef[1:3] and z$coef[4] respectively ## starting value for pfamily in z$coef[5] gnlmm3(y, dist="Student", mu=mu, nest=nest, pmu=c(3.69,-1.19, 0.039), pshape=5, pfamily=0, psd=40, points=3) # or equivalently gnlmm3(y, dist="Student", mu=~exp(b0+b1*sex+b2*dose), nest=nest, pmu=c(3.69,-1.19, 0.039), pshape=5, pfamily=0, psd=40, points=3, envir=reps) ## Not run: # or with identity link print(z <- gnlm::gnlr3(y, dist="Student", mu=~sex+dose, pmu=c(0.1,0,0), pshape=1, pfamily=1)) gnlmm3(y, dist="Student", mu=~sex+dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], pfamily=z$coef[5], psd=50, points=3) # or gnlmm3(y, dist="Student", mu=~b0+b1*sex+b2*dose, nest=nest, pmu=z$coef[1:3], pshape=z$coef[4], pfamily=z$coef[5], psd=50, points=3, envir=reps) # # nonlinear regression with Student t distribution mu <- function(p) p[1]+exp(p[2]+p[3]*sex+p[4]*dose) print(z <- gnlm::gnlr3(y, dist="Student", mu=mu, pmu=c(1,1,0,0), pshape=1, pfamily=1)) gnlmm3(y, dist="Student", mu=mu, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3) # or mu2 <- function(p, linear) p[1]+exp(linear) gnlmm3(y, dist="Student", mu=mu2, linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3) # or gnlmm3(y, dist="Student", mu=~a+exp(linear), linear=~sex+dose, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3) # or gnlmm3(y, dist="Student", mu=~b4+exp(b0+b1*sex+b2*dose), nest=nest, pmu=z$coef[1:4], pshape=z$coef[5], pfamily=z$coef[6], psd=50, points=3, envir=reps) # # include regression for the shape parameter with same mu function shape <- function(p) p[1]+p[2]*sex print(z <- gnlm::gnlr3(y, dist="Student", mu=mu, shape=shape, pmu=z$coef[1:4], pshape=c(z$coef[5],0), pfamily=z$coef[6])) gnlmm3(y, dist="Student", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], pfamily=z$coef[7], psd=5, points=3) # or gnlmm3(y, dist="Student", mu=mu, shape=shape, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], pfamily=z$coef[7], psd=5, points=3, envir=reps) # or gnlmm3(y, dist="Student", mu=~b4+exp(b0+b1*sex+b2*dose), shape=~a1+a2*sex, nest=nest, pmu=z$coef[1:4], pshape=z$coef[5:6], pfamily=z$coef[7], psd=5, points=3, envir=reps) ## End(Not run)
hnlmix fits user-specified nonlinear regression equations to one or
both parameters of the common one and two parameter distributions. One
parameter of the location regression is random with some specified mixing
distribution.
hnlmix( y = NULL, distribution = "normal", mixture = "normal", random = NULL, nest = NULL, mu = NULL, shape = NULL, linear = NULL, pmu = NULL, pshape = NULL, pmix = NULL, prandom = NULL, delta = 1, common = FALSE, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1, eps = 1e-04, points = 5 )hnlmix( y = NULL, distribution = "normal", mixture = "normal", random = NULL, nest = NULL, mu = NULL, shape = NULL, linear = NULL, pmu = NULL, pshape = NULL, pmix = NULL, prandom = NULL, delta = 1, common = FALSE, envir = parent.frame(), print.level = 0, typsize = abs(p), ndigit = 10, gradtol = 1e-05, stepmax = 10 * sqrt(p %*% p), steptol = 1e-05, iterlim = 100, fscale = 1, eps = 1e-04, points = 5 )
y |
A response vector of uncensored data, a two column matrix for
binomial data or censored data, with the second column being the censoring
indicator (1: uncensored, 0: right censored, -1: left censored), or an
object of class, |
distribution |
The distribution for the response: binomial, beta binomial, double binomial, mult(iplicative) binomial, Poisson, negative binomial, double Poisson, mult(iplicative) Poisson, gamma count, Consul generalized Poisson, logarithmic series, geometric, normal, inverse Gauss, logistic, exponential, gamma, Weibull, extreme value, Cauchy, Pareto, Laplace, Levy, beta, simplex, or two-sided power. (For definitions of distributions, see the corresponding [dpqr]distribution help.) |
mixture |
The mixing distribution for the random parameter (whose
initial values are supplied in |
random |
The name of the random parameter in the |
nest |
The cluster variable classifying observations by the unit upon
which they were observed. Ignored if |
mu |
A user-specified formula containing named unknown parameters,
giving the regression equation for the location parameter. This may contain
the keyword, |
shape |
A user-specified formula containing named unknown parameters,
giving the regression equation for the shape parameter. This may contain
the keyword, |
linear |
A formula beginning with ~ in W&R notation, specifying the linear part of the regression function for the location parameter or list of two such expressions for the location and/or shape parameters. |
pmu |
Vector of initial estimates for the location parameters. These must be supplied either in their order of appearance in the formula or in a named list. |
pshape |
Vector of initial estimates for the shape parameters. These must be supplied either in their order of appearance in the expression or in a named list. |
pmix |
If NULL, this parameter is estimated from the variances. If a value is given, it is taken as fixed. |
prandom |
Either one estimate of the random effects or one for each
cluster (see |
delta |
Scalar or vector giving the unit of measurement (always one
for discrete data) for each response value, set to unity by default. For
example, if a response is measured to two decimals, |
common |
If TRUE, the formulae with unknowns for the location and
shape have names in common. All parameter estimates must be supplied in
|
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
typsize |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
steptol |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
fscale |
Arguments for nlm. |
eps |
Arguments for nlm. |
points |
Arguments for nlm. |
It is recommended that initial estimates for pmu and pshape
be obtained from gnlr.
These nonlinear regression models must be supplied as formulae where
parameters are unknowns. (See finterp.)
A list of class hnlmix is returned that contains all of the
relevant information calculated, including error codes.
The two variances and shrinkage estimates of the random effects are provided.
J.K. Lindsey
dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8) #y <- rgamma(20,2+0.3*dose,scale=2)+rep(rnorm(4,0,4),rep(5,4)) y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699, 4.359469, 1.900681, 17.425948, 4.503345, 2.691792, 5.731100, 10.534971, 11.220260, 6.968932, 4.094357, 16.393806, 14.656584, 8.786133, 20.972267, 17.178012) resp <- restovec(matrix(y, nrow=4, byrow=TRUE), name="y") reps <- rmna(resp, tvcov=tvctomat(matrix(dose, nrow=4, byrow=TRUE), name="dose")) # same linear normal model with random normal intercept fitted four ways # compare with growth::elliptic(reps, model=~dose, preg=c(0,0.6), pre=4) glmm(y~dose, nest=individuals, data=reps) gnlmm(reps, mu=~dose, pmu=c(8.7,0.25), psh=3.5, psd=3) hnlmix(reps, mu=~a+b*dose+rand, random="rand", pmu=c(8.7,0.25), pshape=3.44, prandom=0) # gamma model with log link and random normal intercept fitted three ways glmm(y~dose, family=Gamma(link=log), nest=individuals, data=reps, points=8) gnlmm(reps, distribution="gamma", mu=~exp(a+b*dose), pmu=c(2,0.03), psh=1, psd=0.3) hnlmix(reps, distribution="gamma", mu=~exp(a+b*dose+rand), random="rand", pmu=c(2,0.04), pshape=1, prandom=0) # gamma model with log link and random gamma mixtures hnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a*rand+b*dose), random="rand", pmu=c(2,0.04), pshape=1.24, prandom=1) hnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a+b*dose)*rand, random="rand", pmu=c(2,0.04), pshape=1.24, prandom=1)dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8) #y <- rgamma(20,2+0.3*dose,scale=2)+rep(rnorm(4,0,4),rep(5,4)) y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699, 4.359469, 1.900681, 17.425948, 4.503345, 2.691792, 5.731100, 10.534971, 11.220260, 6.968932, 4.094357, 16.393806, 14.656584, 8.786133, 20.972267, 17.178012) resp <- restovec(matrix(y, nrow=4, byrow=TRUE), name="y") reps <- rmna(resp, tvcov=tvctomat(matrix(dose, nrow=4, byrow=TRUE), name="dose")) # same linear normal model with random normal intercept fitted four ways # compare with growth::elliptic(reps, model=~dose, preg=c(0,0.6), pre=4) glmm(y~dose, nest=individuals, data=reps) gnlmm(reps, mu=~dose, pmu=c(8.7,0.25), psh=3.5, psd=3) hnlmix(reps, mu=~a+b*dose+rand, random="rand", pmu=c(8.7,0.25), pshape=3.44, prandom=0) # gamma model with log link and random normal intercept fitted three ways glmm(y~dose, family=Gamma(link=log), nest=individuals, data=reps, points=8) gnlmm(reps, distribution="gamma", mu=~exp(a+b*dose), pmu=c(2,0.03), psh=1, psd=0.3) hnlmix(reps, distribution="gamma", mu=~exp(a+b*dose+rand), random="rand", pmu=c(2,0.04), pshape=1, prandom=0) # gamma model with log link and random gamma mixtures hnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a*rand+b*dose), random="rand", pmu=c(2,0.04), pshape=1.24, prandom=1) hnlmix(reps, distribution="gamma", mixture="gamma", mu=~exp(a+b*dose)*rand, random="rand", pmu=c(2,0.04), pshape=1.24, prandom=1)
kalcount is designed to handle repeated measurements models with
time-varying covariates. The distributions have two extra parameters as
compared to the functions specified by intensity and are generally
longer tailed than those distributions. Dependence among observations on a
unit can be through gamma or power variance family frailties (a type of
random effect), with or without autoregression, or serial dependence over
time.
kalcount( response = NULL, times = NULL, origin = 0, intensity = "exponential", depend = "independence", update = "Markov", mu = NULL, shape = NULL, density = FALSE, ccov = NULL, tvcov = NULL, preg = NULL, ptvc = NULL, pbirth = NULL, pintercept = NULL, pshape = NULL, pinitial = 1, pdepend = NULL, pfamily = NULL, envir = parent.frame(), print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )kalcount( response = NULL, times = NULL, origin = 0, intensity = "exponential", depend = "independence", update = "Markov", mu = NULL, shape = NULL, density = FALSE, ccov = NULL, tvcov = NULL, preg = NULL, ptvc = NULL, pbirth = NULL, pintercept = NULL, pshape = NULL, pinitial = 1, pdepend = NULL, pfamily = NULL, envir = parent.frame(), print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )
response |
A list of two column matrices with counts and corresponding
times for each individual, one matrix or dataframe of counts, or an object
of class, |
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,
|
origin |
If the time origin is to be before the start of observations, a positive constant to be added to all times. |
intensity |
The form of function to be put in the Pareto distribution. Choices are exponential, Weibull, gamma, log normal, log logistic, log Cauchy, log Student, and gen(eralized) logistic. |
depend |
Type of dependence. Choices are |
update |
Type of for serial dependence. Choices are |
mu |
A regression function for the location parameter or a formula
beginning with ~, specifying either a linear regression function in the
Wilkinson and Rogers notation (a log link is assumed) or a general function
with named unknown parameters. Give the initial estimates in |
shape |
A regression function for the shape parameter or a formula beginning with ~, specifying either a linear regression function in the Wilkinson and Rogers notation or a general function with named unknown parameters. It must yield one value per observation. |
density |
If TRUE, the density of the function specified in
|
ccov |
A vector or matrix containing time-constant baseline covariates
with one row per individual, a model formula using vectors of the same
size, or an object of class, |
tvcov |
A list of matrices with time-varying covariate values,
observed in the time periods in |
preg |
Initial parameter estimates for the regression model: intercept
plus one for each covariate in |
ptvc |
Initial parameter estimates for the coefficients of the
time-varying covariates, as many as in |
pbirth |
If supplied, this is the initial estimate for the coefficient of the birth model. |
pintercept |
The initial estimate of the intercept for the generalized logistic intensity. |
pshape |
An initial estimate for the shape parameter of the intensity
function (except exponential intensity). If |
pinitial |
An initial estimate for the initial parameter. With
|
pdepend |
An initial estimate for the serial dependence parameter. For
|
pfamily |
An optional initial estimate for the second parameter of a
two-parameter power variance family mixture instead of the default gamma
mixture. This yields a gamma mixture as |
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
steptol |
Arguments for nlm. |
fscale |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
typsize |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
By default, a gamma mixture of the distribution specified in
intensity is used, as the conditional distribution in the
serial dependence models, and as a symmetric multivariate (random
effect) model for frailty dependence.
Unless specified otherwise, the time origin is taken to be zero. The given times are the ends of the periods in which the counts occurred.
Here, the variance, with exponential intensity, is a quadratic function of
the mean, whereas, for nbkal, it is proportional to
the mean function.
If the counts on a unit are clustered, not longitudinal, use the failty dependence with the default exponential intensity, yielding a multivariate negative binomial distribution.
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 finterp.)
Marginal and individual profiles can be plotted using
mprofile and iprofile and
residuals with plot.residuals.
A list of classes kalcount and recursive is returned.
J.K. Lindsey
treat <- c(0,0,1,1) tr <- tcctomat(treat) dose <- matrix(c(9,13,16,7,12,6,9,10,11,9,10,10,7,9,9,9,8,10,15,4), ncol=5,byrow=TRUE) dd <- tvctomat(dose) y <- restovec(structure(c(6, 4, 0, 0, 3, 6, 1, 1, 1, 5, 0, 0, 0, 4, 0, 1, 0, 13, 0, 3), dim = 4:5)) reps <- rmna(y, ccov=tr, tvcov=dd) kalcount(y, intensity="log normal", dep="independence", preg=0.3, pshape=0.1) kalcount(y, intensity="log normal", dep="frailty", pinitial=0.1, preg=1, psh=0.1) kalcount(y, intensity="log normal", dep="serial", pinitial=0.1, preg=1, pdep=0.75, psh=0.1) # random effect and autoregression (commented out: AR difficult to estimate) #kalcount(y, intensity="log normal", dep="frailty", pinitial=0.1, # pdep=0.5, preg=1, psh=0.1) # add time-constant variable kalcount(y, intensity="log normal", pinitial=1, psh=1, preg=c(0.8,0.11), ccov=treat) # or kalcount(y, intensity="log normal", mu=~b0+b1*treat, pinitial=1, psh=.1, preg=c(0.4,-0.04), envir=reps) # add time-varying variable kalcount(y, intensity="log normal", pinitial=1, psh=1, preg=c(-1,2), ccov=treat, ptvc=0, tvc=dose) # or equivalently, from the environment dosev <- as.vector(t(dose)) kalcount(y, intensity="log normal", mu=~b0+b1*rep(treat,rep(5,4))+b2*dosev, pinitial=1, psh=1, ptvc=c(-1,2,0)) # or from the reps data object kalcount(y, intensity="log normal", mu=~b0+b1*treat+b2*dose, pinitial=1, psh=1, ptvc=c(-1,2,0), envir=reps) # try power variance family kalcount(y, intensity="log normal", mu=~b0+b1*treat+b2*dose, pinitial=1, psh=1, ptvc=c(-1,2,0.1), envir=reps, pfamily=0.8)treat <- c(0,0,1,1) tr <- tcctomat(treat) dose <- matrix(c(9,13,16,7,12,6,9,10,11,9,10,10,7,9,9,9,8,10,15,4), ncol=5,byrow=TRUE) dd <- tvctomat(dose) y <- restovec(structure(c(6, 4, 0, 0, 3, 6, 1, 1, 1, 5, 0, 0, 0, 4, 0, 1, 0, 13, 0, 3), dim = 4:5)) reps <- rmna(y, ccov=tr, tvcov=dd) kalcount(y, intensity="log normal", dep="independence", preg=0.3, pshape=0.1) kalcount(y, intensity="log normal", dep="frailty", pinitial=0.1, preg=1, psh=0.1) kalcount(y, intensity="log normal", dep="serial", pinitial=0.1, preg=1, pdep=0.75, psh=0.1) # random effect and autoregression (commented out: AR difficult to estimate) #kalcount(y, intensity="log normal", dep="frailty", pinitial=0.1, # pdep=0.5, preg=1, psh=0.1) # add time-constant variable kalcount(y, intensity="log normal", pinitial=1, psh=1, preg=c(0.8,0.11), ccov=treat) # or kalcount(y, intensity="log normal", mu=~b0+b1*treat, pinitial=1, psh=.1, preg=c(0.4,-0.04), envir=reps) # add time-varying variable kalcount(y, intensity="log normal", pinitial=1, psh=1, preg=c(-1,2), ccov=treat, ptvc=0, tvc=dose) # or equivalently, from the environment dosev <- as.vector(t(dose)) kalcount(y, intensity="log normal", mu=~b0+b1*rep(treat,rep(5,4))+b2*dosev, pinitial=1, psh=1, ptvc=c(-1,2,0)) # or from the reps data object kalcount(y, intensity="log normal", mu=~b0+b1*treat+b2*dose, pinitial=1, psh=1, ptvc=c(-1,2,0), envir=reps) # try power variance family kalcount(y, intensity="log normal", mu=~b0+b1*treat+b2*dose, pinitial=1, psh=1, ptvc=c(-1,2,0.1), envir=reps, pfamily=0.8)
kalseries is designed to handle repeated measurements models with
time-varying covariates. The distributions have two extra parameters as
compared to the functions specified by intensity and are generally
longer tailed than those distributions. Dependence among observations on a
unit can be through gamma or power variance family frailties (a type of
random effect), with or without autoregression, or one of two types of
serial dependence over time.
kalseries( response = NULL, times = NULL, intensity = "exponential", depend = "independence", mu = NULL, shape = NULL, density = FALSE, ccov = NULL, tvcov = NULL, torder = 0, interaction = NULL, preg = NULL, ptvc = NULL, pintercept = NULL, pshape = NULL, pinitial = 1, pdepend = NULL, pfamily = NULL, delta = NULL, transform = "identity", link = "identity", envir = parent.frame(), print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )kalseries( response = NULL, times = NULL, intensity = "exponential", depend = "independence", mu = NULL, shape = NULL, density = FALSE, ccov = NULL, tvcov = NULL, torder = 0, interaction = NULL, preg = NULL, ptvc = NULL, pintercept = NULL, pshape = NULL, pinitial = 1, pdepend = NULL, pfamily = NULL, delta = NULL, transform = "identity", link = "identity", envir = parent.frame(), print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )
response |
A list of two column matrices with responses and
corresponding times for each individual, one matrix or dataframe of
response values, or an object of class, |
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,
|
intensity |
The form of function to be put in the Pareto distribution. Choices are exponential, Weibull, gamma, normal, logistic, Cauchy, log normal, log logistic, log Cauchy, log Student, inverse Gauss, and gen(eralized) logistic. (For definitions of distributions, see the corresponding [dpqr]distribution help.) |
depend |
Type of dependence. Choices are |
mu |
A regression function for the location parameter or a formula
beginning with ~, specifying either a linear regression function in the
Wilkinson and Rogers notation or a general function with named unknown
parameters. Give the initial estimates in |
shape |
A regression function for the shape parameter or a formula beginning with ~, specifying either a linear regression function in the Wilkinson and Rogers notation or a general function with named unknown parameters. It must yield one value per observation. |
density |
If TRUE, the density of the function specified in
|
ccov |
A vector or matrix containing time-constant baseline covariates
with one row per individual, a model formula using vectors of the same
size, or an object of class, |
tvcov |
A list of matrices with time-varying covariate values,
observed at the event times in |
torder |
The order of the polynomial in time to be fitted. |
interaction |
Vector of length equal to the number of time-constant
covariates, giving the levels of interactions between them and the
polynomial in time in the |
preg |
Initial parameter estimates for the regression model:
intercept, one for each covariate in |
ptvc |
Initial parameter estimates for the coefficients of the
time-varying covariates, as many as in |
pintercept |
The initial estimate of the intercept for the generalized logistic intensity. |
pshape |
An initial estimate for the shape parameter of the intensity
function (except exponential intensity). If |
pinitial |
An initial estimate for the initial parameter. With
|
pdepend |
An initial estimate for the serial dependence parameter. For
|
pfamily |
An optional initial estimate for the second parameter of a
two-parameter power variance family mixture instead of the default gamma
mixture. This yields a gamma mixture as |
delta |
Scalar or vector giving the unit of measurement for each
response value, set to unity by default. For example, if a response is
measured to two decimals, delta=0.01. 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, |
transform |
Transformation of the response variable: |
link |
Link function for the mean: |
envir |
Environment in which model formulae are to be interpreted or a
data object of class, |
print.level |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
steptol |
Arguments for nlm. |
fscale |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
typsize |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
By default, a gamma mixture of the distribution specified in
intensity is used, as the conditional distribution in the
Markov and serial dependence models, and as a symmetric
multivariate (random effect) model for frailty dependence. For
example, with a Weibull intensity and frailty dependence,
this yields a multivariate Burr distribution and with Markov or
serial dependence, univariate Burr conditional distributions.
If a value for pfamily is used, the gamma mixture is replaced by a
power variance family mixture.
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 finterp.)
Marginal and individual profiles can be plotted using
mprofile and iprofile and
residuals with plot.residuals.
A list of classes kalseries and recursive is
returned.
J.K. Lindsey
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) ## Not run: # 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)) ## End(Not run)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) ## Not run: # 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)) ## End(Not run)
marg.hom fits a marginal homogeneity model to a contingency table
that has two margins of equal size.
marg.hom(freq, marg1, marg2)marg.hom(freq, marg1, marg2)
freq |
Vector of frequencies |
marg1 |
Factor variable for the first margin |
marg2 |
Factor variable for the second margin |
A list containing the call, the model, the deviance, the degrees of freedom, the aic, the fitted values, and the residuals is returned.
J.K. Lindsey
# 4x4x2 table in Freq, with margins indexed by Left and Right Freq <- rpois(32,10) Left <- gl(4,1,32) Right <- gl(4,4,32) marg.hom(Freq, Left, Right)# 4x4x2 table in Freq, with margins indexed by Left and Right Freq <- rpois(32,10) Left <- gl(4,1,32) Right <- gl(4,4,32) marg.hom(Freq, Left, Right)
nbkal fits a negative binomial regression with Kalman update over
time. The variance is proportional to the mean function, whereas, for
kalcount with exponential intensity, it is a
quadratic function of the mean.
nbkal( response, times, mu, preg, pdepend, kalman = TRUE, print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )nbkal( response, times, mu, preg, pdepend, kalman = TRUE, print.level = 0, ndigit = 10, gradtol = 1e-05, steptol = 1e-05, fscale = 1, iterlim = 100, typsize = abs(p), stepmax = 10 * sqrt(p %*% p) )
response |
A list of two column matrices with counts and corresponding
times for each individual, one matrix or dataframe of counts, or an object
of class, response (created by |
times |
When response is a matrix, a vector of possibly unequally spaced times when they are the same for all individuals or a matrix of times. Not necessary if equally spaced. Ignored if response has class, response or repeated. |
mu |
The mean function. |
preg |
The initial parameter estimates for the mean function. |
pdepend |
The estimates for the dependence parameters, either one or three. |
kalman |
If TRUE, fits the kalman update model, otherwise, a standard negative binomial distribution. |
print.level |
Arguments for nlm. |
ndigit |
Arguments for nlm. |
gradtol |
Arguments for nlm. |
steptol |
Arguments for nlm. |
fscale |
Arguments for nlm. |
iterlim |
Arguments for nlm. |
typsize |
Arguments for nlm. |
stepmax |
Arguments for nlm. |
Marginal and individual profiles can be plotted using
mprofile and iprofile and
residuals with plot.residuals.
A list of classes nbkal and recursive is returned.
P. Lambert and J.K. Lindsey
Lambert, P. (1996) Applied Statistics 45, 31-38.
Lambert, P. (1996) Biometrics 52, 50-55.
y <- matrix(rnbinom(20,5,0.5), ncol=5) times <- matrix(rep(seq(10,50,by=10),4), ncol=5, byrow=TRUE) y0 <- matrix(rep(rnbinom(5,5,0.5),4), ncol=5, byrow=TRUE) mu <- function(p) p[1]*log(y0)+(times<30)*p[2]* (times-30)+(times>30)*p[3]*(times-30) nbkal(y, preg=c(1.3,0.008,-0.05), times=times, pdep=1.2, mu=mu)y <- matrix(rnbinom(20,5,0.5), ncol=5) times <- matrix(rep(seq(10,50,by=10),4), ncol=5, byrow=TRUE) y0 <- matrix(rep(rnbinom(5,5,0.5),4), ncol=5, byrow=TRUE) mu <- function(p) p[1]*log(y0)+(times<30)*p[2]* (times-30)+(times>30)*p[3]*(times-30) nbkal(y, preg=c(1.3,0.008,-0.05), times=times, pdep=1.2, mu=mu)