Segmented relationships in linear mixed models

Fits (generalized) linear mixed models with a segmented relationship between the response and a numeric covariate. Random effects are allowed in each model parameter, including the breakpoint.

Usage

# S3 method for class 'lme'
segmented(obj, seg.Z, psi, npsi = 1, fixed.psi = NULL, 
    control = seg.control(), model = TRUE, 
    z.psi = ~1, x.diff = ~1, random = NULL, 
    random.noG = NULL, start.pd = NULL, psi.link = c("identity", "logit"), 
    start = NULL, data, fixed.parms = NULL,...)

Arguments

obj

A fit returned by lme() or simply its call (if using segmented.lme). See example below. This represents the linear mixed model where the segmented relationship is added. It has to include the argument data.

seg.Z

A one-sided formula indicating the segmented variable, i.e. the quantitative variable having a segmented relationship with the response. In longitudinal studies typically it is the time.

psi

An optional starting value for the breakpoint. If missing, a starting value is obtained via the nadir estimate of a quadratic fit. When provided it may be a single numeric value or a vector of length equal to the number of clusters (i.e. subjects).

z.psi

Optional. A one-sided formula meaning the covariates in the sub-regression model for the changepoint parameter. Default to ~1.

x.diff

Optional. A one-sided formula meaning the covariates in the sub-regression model for the difference-in-slopes parameter. Default to ~1 for no covariate for the difference-in-slope parameter.

npsi

Ignored. Currently only npsi=1 is allowed.

fixed.psi

Ignored.

control

A list returned by seg.control, in particular display, n.boot for the bootstrap restarting.

model

Ignored.

random

A list, as the one supplied in random of lme() including the random effects. Default to NULL, meaning that the same random effect structure of the initial lme fit supplied in obj should be used. When specified, this list could include the variables 'G0' and 'U'. G0 means random effects in the breakpoints and U means random effects in the slope-difference parameter. Assuming id is the the cluster variable and x the segmented variable, some examples are

random = list(id = pdDiag(~1 + x + U)) #ind. random eff. (changepoint fixed)

random = list(id = pdDiag(~1 + x + U + G0)) #ind. random eff. (in the changepoint too)

random = list(id=pdBlocked(list(pdSymm(~1+x), pdSymm(~U+G0-1)))) #block diagonal

random.noG

Ignored.

start.pd

An optional starting value for the variances of the random effects. It should be coherent with the specification in random.

psi.link

The link function used to specify the sub-regression model for the breakpoint $psi$. The identity (default) assumes $$\psi_i=\eta_i$$ while the logit link is $$\psi_i=(m+M*exp(\eta_i))/(1+exp(\eta_i))$$ where $m$ and $M$ are the observed minimum and maximum of the segmented variable in seg.Z. In each case the `linear predictor' is $\eta_i=\kappa_0+z_i^T\kappa_1+k_i$, where $z^T$ includes the covariates specified in z.psi and the $k_i$s are the changepoint random effects included by means of G0 in the random argument.

start

An optional list including the starting values for the difference-in-slopes parameter, delta0 and delta, and the changepoint parameter, kappa and kappa0. When provided, 'kappa0' overwrites 'psi'.

If provided, the components 'delta' and 'kappa' should be named vectors with length and names matching length and names in x.diff and z.psi respectively. The component delta0 can be a scalar or a vector with length equal to the number of clusters (subjects).

data

the dataframe where the variables are stored. If missing, the dataframe of the "lme" (or "glmmPQL" ) fit obj is assumed.

fixed.parms

An optional named vector representing the coefficients of the changepoint to be maintained fixed during the estimation process. Allowed names are "G0" or any variable (in the dataframe) supposed to affect the location of breakpoints. For instance fixed.parms=c(G0=.3) implies a fixed value for the changepoint. Notice if you use the same variable in fixed.parms and in z.psi, for instance fixed.parms=c(x2=.3) and z.psi=~x2, a warning is printed and the coefficient "G.x2" is estimated to maximize the log likelihood given that fixed value. As an example, suppose the unconstrained estimated coefficient for x2, say, in z.psi is 0.5; if in a new call both fixed.parms=c(x2=.4) and z.psi=~x2 are included, the estimate of "G.x2" will be (approximately) 0.1. Essentially, if you really want to fix the parameters in fixed.parms, then do not include the same covariates in z.psi.

...

Ignored

Details

The function fits segmented mixed regression models, i.e. segmented models with random effects also in the slope-difference and change-point parameters.

Value

A list of class segmented.lme with several components. The most relevant are

lme.fit: The fitted lme object at convergence
lme.fit.noG: The fitted lme object at convergence assuming known the breakpoints
psi.i: The subject/cluster-specific change points (fixed + random). It includes 2 attributes: attr(,"ni") for the number of measurements in each 'cluster', and attr(,"is.break") a vector of logicals indicating if the breakpoint for each subject i can be reliable (TRUE) or not (FALSE). Here 'reliable' simply means within the covariate range (for subject i). See also argument nq.
fixed.eta.psi: The fixed-effect linear predictor for the change points regression equation. These values will different among 'clusters' only if at least one covariate has been specified in z.psi.
fixed.eta.delta: The fixed-effect linear predictor of the slope difference regression equation. These values will different among 'clusters' only if at least one covariate has been specified in x.diff.

References

Muggeo V., Atkins D.C., Gallop R.J., Dimidjian S. (2014) Segmented mixed models with random changepoints: a maximum likelihood approach with application to treatment for depression study. Statistical Modelling, 14, 293-313.

Muggeo V. (2016) Segmented mixed models with random changepoints in R. Working paper available on RG. doi: 10.13140/RG.2.1.4180.8402

Author

Vito M.R. Muggeo vito.muggeo@unipa.it

Note

segmented.lme() relies on lme() using REML by default.

If fit is the segmented.lme fit, use VarCorr(fit$lme.fit) to extract the random effect covariance matrix.

Warning

The function deals with estimation of a single breakpoint only (with or without random effects).

Examples


if (FALSE) { # \dontrun{
library(nlme)
data(Cefamandole)
Cefamandole$lTime <-log(Cefamandole$Time)
Cefamandole$lconc <-log(Cefamandole$conc)

o<-lme(lconc ~ lTime, random=~1|Subject, data=Cefamandole)

os<-segmented.lme(o, ~lTime, random=list(Subject=pdDiag(~1+lTime+U+G0)), 
  control=seg.control(n.boot=0, display=TRUE))
slope(os)


####################################################
# covariate effect on the changepoint and slope diff


#let's assume a new subject-specific covariates..
set.seed(69)
Cefamandole$z <- rep(runif(6), rep(14,6))
Cefamandole$group <- gl(2,42,labels=c('a','b'))

#Here 'group' affects the slopes and 'z' affects the changepoint 

o1 <-lme(lconc ~ lTime*group, random=~1|Subject, data=Cefamandole)
os1 <- segmented(o1, ~lTime, x.diff=~group, z.psi=~z, 
  random=list(Subject=pdDiag(~1+lTime+U+G0)))

slope(os1, by=list(group="a")) #the slope estimates in group="a" (baseline level)
slope(os1, by=list(group="b")) #the slope estimates in group="b" 


###################################################
# A somewhat "complicated" example:
#     i)  strong heterogeneity in the changepoints
#     ii) No changepoint for the Subject #7 (added) 

d<-Cefamandole
d$x<- d$lTime
d$x[d$Subject==1]<- d$lTime[d$Subject==1]+3
d$x[d$Subject==5]<- d$lTime[d$Subject==5]+5
d$x[d$Subject==3]<- d$lTime[d$Subject==3]-5
d<-rbind(d, d[71:76,])
d$Subject <- factor(d$Subject, levels=c(levels(d$Subject),"7")) 
d$Subject[85:90] <- rep("7",6)

o<-lme(lconc ~ x, random=~1|Subject, data=d)
os2<-segmented.lme(o, ~x, random=list(Subject=pdDiag(~1+x+U+G0)), 
  control=seg.control(n.boot=5, display=TRUE))

#plots with common x- and y- scales (to note heterogeneity in the changepoints)
plot(os2, n.plot = c(3,3)) 
os2$psi.i
attr(os2$psi.i, "is.break") #it is FALSE for Subject #7

#plots with subject-specific scales
plot(os2, n.plot = c(3,3), xscale=-1, yscale = -1) 
} # }

Segmented relationships in linear mixed models

Usage

Arguments

Details

Value

References

Author

Note

Warning

See also

Examples