The brokenstick()
function fits an irregularly observed series
of measurements onto a userspecified grid of points.
The model codes de grid by a series of linear Bsplines.
Differences between observations are expressed by one random
effect per grid point. When multiple set of series are modelled,
each modeled trajectory consists of straight lines that join at the
chosen grid points, and hence look like a broken stick.
brokenstick(x, ...) # S3 method for default brokenstick(x, ...) # S3 method for formula brokenstick( formula, data, ..., knots = NULL, boundary = NULL, k = NULL, degree = 1L, method = c("lmer", "kr"), control = control_brokenstick(), seed = NA ) # S3 method for data.frame brokenstick( x, y, group, ..., knots = NULL, boundary = NULL, k = NULL, degree = 1L, method = c("lmer", "kr"), control = control_brokenstick(), seed = NA ) # S3 method for matrix brokenstick( x, y, group, ..., knots = NULL, boundary = NULL, k = NULL, degree = 1L, method = c("lmer", "kr"), control = control_brokenstick(), seed = NA ) # S3 method for numeric brokenstick( x, y, group, ..., knots = NULL, boundary = NULL, k = NULL, degree = 1L, method = c("lmer", "kr"), control = control_brokenstick(), seed = NA )
x  Predictor variables. Depending on the context:
If 

...  Not currently used, but required for extensibility. 
formula  A formula specifying the outcome terms on the
lefthand side, the predictor term on the righthand side and
the group variable after the Note: This formula specification is specific to the 
data  When a formula is used,

knots  Optional, but recommended. Numerical vector with the
locations of the breaks to be placed on the values of the predictor.
Values outside the range of the data will extend the 
boundary  Optional. Numerical vector of
length 2 with the left and right boundary knots. The 
k  Optional, a convenience parameter for the number of
internal knots. If specified, then Note: Knots specified via 
degree  the degree of the spline. The broken stick model
requires linear splines, so the default is 
method  Estimation method. Either 
control  A list with parameters. Use 
seed  Seed number for 
y  Outcome variable. When

group  Grouping variable. When

A brokenstick
object.
The variancecovariance matrix of the random effects absorbs the
relations over time. By default, this matrix is estimated
as unstructured by lme4::lmer()
.
This estimate may be unstable if
the number of children is small relative to the number of specified
knots. The function can be time consuming for data sets with thousands of
children.
An alternative  often faster for models with many random effects 
is to use the Bayesian KasimRaudenbush sampler (method kr
). That
method also allow for enforcing a simple structure on the
variancecovariance matrix of the random effects. Currently, there
are three such models: argyle
, cole
and none
. See kr()
and
control_brokenstick()
for more details.
train < smocc_200[1:1198, ] # fit with implicit boundary c(0, 3) fit < brokenstick(hgt.z ~ age  id, data = train, knots = 0:3) # using KR sampler fit < brokenstick(hgt.z ~ age  id, data = train, knots = 0:3, method = "kr", seed = 1) # \donttest{ knots < round(c(0, 1, 2, 3, 6, 9, 12, 15, 18, 24) / 12, 4) boundary < c(0, 3) fit_lmer < brokenstick(hgt.z ~ age  id, data = train, knots = knots, boundary = boundary)#> Warning: number of observations (=1185) <= number of random effects (=1364) for term (0 + age_0 + age_0.0833 + age_0.1667 + age_0.25 + age_0.5 + age_0.75 + age_1 + age_1.25 + age_1.5 + age_2 + age_3  id); the randomeffects parameters and the residual variance (or scale parameter) are probably unidentifiable#>fit_kr < brokenstick(hgt.z ~ age  id, data = train, knots = knots, boundary = boundary, method = "kr") # } # Four ways to specify the same model # Formula interface mod1 < brokenstick(hgt.z ~ age  id, train) # XY interface  numeric vector mod2 < with(train, brokenstick(age, hgt.z, id)) identical(mod1, mod2)#> [1] FALSE# XY interface  data.frame mod3 < with(train, brokenstick(data.frame(age), hgt.z, id)) identical(mod1, mod3)#> [1] FALSE# XY interface  matrix tt < as.matrix(train[, c(1, 2, 7)]) mod4 < brokenstick(tt[, "age", drop = FALSE], tt[, "hgt.z", drop = FALSE], tt[, "id", drop = FALSE]) identical(mod1, mod4)#> [1] FALSE