Skip to content

Commit 3e8da94

Browse files
Merge pull request #6 from growthcharts/revision
brokenstick 2.0.0
2 parents 9b969af + 971f21a commit 3e8da94

File tree

161 files changed

+4900
-3261
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

161 files changed

+4900
-3261
lines changed

.gitignore

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,4 +27,6 @@ vignettes/*.R
2727
# Temporary files created by R markdown
2828
*.utf8.md
2929
*.knit.md
30-
.Rproj.user
30+
.Rproj.user
31+
/doc/
32+
/Meta/

DESCRIPTION

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ Description: The broken stick model describes a set of individual curves by a
1515
Depends:
1616
R (>= 3.5.0)
1717
Imports:
18+
coda,
1819
dplyr,
1920
lme4,
2021
matrixsampling,
@@ -45,3 +46,4 @@ License: MIT + file LICENSE
4546
LazyData: TRUE
4647
VignetteBuilder: knitr
4748
Roxygen: list(markdown = TRUE)
49+
RoxygenNote: 7.1.2

NAMESPACE

Lines changed: 9 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,24 @@
11
# Generated by roxygen2: do not edit by hand
22

3-
S3method(brokenstick,data.frame)
4-
S3method(brokenstick,default)
5-
S3method(brokenstick,formula)
6-
S3method(brokenstick,matrix)
7-
S3method(brokenstick,numeric)
3+
S3method(coef,brokenstick)
84
S3method(fitted,brokenstick)
95
S3method(model.frame,brokenstick)
106
S3method(model.matrix,brokenstick)
117
S3method(plot,brokenstick)
128
S3method(predict,brokenstick)
139
S3method(print,brokenstick)
10+
S3method(print,summary.brokenstick)
1411
S3method(residuals,brokenstick)
15-
export(EB)
12+
S3method(summary,brokenstick)
1613
export(brokenstick)
17-
export(control_brokenstick)
1814
export(control_kr)
1915
export(get_knots)
16+
export(get_omega)
2017
export(get_r2)
2118
export(kr)
22-
export(make_basis)
23-
export(parse_formula)
19+
export(plot_trajectory)
20+
export(set_control)
21+
importFrom(coda,mcmc)
2422
importFrom(dplyr,`%>%`)
2523
importFrom(dplyr,bind_cols)
2624
importFrom(dplyr,bind_rows)
@@ -36,14 +34,15 @@ importFrom(lme4,VarCorr)
3634
importFrom(lme4,fixef)
3735
importFrom(lme4,lmer)
3836
importFrom(lme4,lmerControl)
39-
importFrom(lme4,ranef)
37+
importFrom(lme4,ngrps)
4038
importFrom(matrixsampling,rwishart)
4139
importFrom(methods,slot)
4240
importFrom(rlang,.data)
4341
importFrom(rlang,arg_match)
4442
importFrom(splines,bs)
4543
importFrom(stats,approx)
4644
importFrom(stats,as.formula)
45+
importFrom(stats,coef)
4746
importFrom(stats,cor)
4847
importFrom(stats,cov2cor)
4948
importFrom(stats,fitted)

NEWS.md

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,35 @@
1+
# brokenstick 2.0.0
2+
3+
## Main changes
4+
5+
1. Function `brokenstick()` in version `2.0.0` sets the Kasim-Raudenbush sampler as the default method. The former method `lme4::lmer()` remains available by setting `method = "lmer"` argument.
6+
7+
2. Version `2.0.0` adopts the variable names of the `coda` package (e.g., `start`, `end`, `thin`, `niter`, and so on) and stores the results of the Kasim-Raudenbush sampler as objects of class `mcmc`.
8+
9+
3. For `method = "kr"` one may now inspect the solution of the sampler by standard functions from the `coda` package. For `method = "lmer"` we can apply functions from the `lme4` package for `merMod` objects.
10+
11+
4. Version `2.0.0` redefines the `brokenstick` class. New entries include `call`, `formula`, `internal`, `sample`, `light`, `data`, `imp` and `mod`. Removed entries are `knots` (renamed to `internal`) and `draws` (renamed to `imp`). We may omit the `newdata` argument for the training data. Setting `light = TRUE` creates a small version of the `brokenstick` object. Objects of class `brokenstick` are not backwards compatible, so one should regenerate objects of class `brokenstick` in order use newer features in `2.0.0`.
12+
13+
5. Version `2.0.0` conforms to classic model fitting interface in `R`. Renames the `new_data` argument to `newdata` to conform to `predict.lm()`. Methods `plot()` and `predict()` no longer require a `newdata` argument. All special cases of `predict()` updated and explained in documentation and examples.
14+
15+
6. Version `2.0.0` adds methods `coef()`, `fitted()`, `model.frame()`, `model.matrix()`, `print()` and `summary` for the `brokenstick` object.
16+
17+
7. Simplifies algorithmic control. Renames `control_brokenstick()` to `set_control()` and removes a layer in the control list.
18+
19+
## Minor changes
20+
21+
- Stabilises the `rgamma()` calls in KR-algorithm for edge cases.
22+
- `predict_brokenstick()` can now work with the both (internal) training and (external) test data.
23+
- Removes the superfluous `type` argument from `predict.brokenstick()`
24+
- Adds a function `get_omega()` to extract the variance-covariance matrix of the broken stick estimates
25+
- Adds choice `"dropfirst"` to `get_knots()`
26+
- Improves error messages of edge cases in `test-brokenstick_edge.R`
27+
- Perform stricter tests on arguments of `brokenstick()`
28+
- Introduces argument `warn_splines` in `make_basis()` to suppress uninteresting warns from `splines::bs()`
29+
- Removes superfluous `knotnames` argument in `make_basis()`
30+
- Argument `x` in `make_basis()` is now a vector instead of a column vector
31+
- Introduces new `xname` argument in `make_basis()` to set the xname
32+
133
# brokenstick 1.1.1
234

335
- Handles an edge case that crashed `predict()`

R/EB.R

Lines changed: 5 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -23,22 +23,10 @@
2323
#' Skrondal, A., Rabe-Hesketh, S. (2009).
2424
#' Prediction in multilevel generalized linear models.
2525
#' J. R. Statist. Soc. A, 172, 3, 659-687.
26-
#' @examples
27-
#' #
28-
#' # EB estimate random effect for child id 10001
29-
#' data <- smocc_200[smocc_200$id == 10001, ]
30-
#' y <- data$hgt.z
31-
#' X <- make_basis(data$age,
32-
#' knots = fit_200$knots,
33-
#' boundary = fit_200$boundary)
34-
#' EB(fit_200, y, X)
35-
#' @export
3626
EB <- function(model, y, X, Z = X, BS = TRUE) {
3727

38-
if (!inherits(model, "brokenstick")) stop("Argument `model` not of class brokenstick.")
39-
40-
# X should be a matrix
41-
if (!is.matrix(X)) stop("Argument 'X' is not a matrix.")
28+
stopifnot(inherits(model, "brokenstick"),
29+
is.matrix(X))
4230

4331
# eliminate missing outcomes
4432
select <- !(is.na(y) | is.na(X[, 1]))
@@ -56,12 +44,10 @@ EB <- function(model, y, X, Z = X, BS = TRUE) {
5644
X <- matrix(X[select, ], ncol = dim(X)[2])
5745
beta <- matrix(model$beta, ncol = 1)
5846

59-
# construct appropriate matrices
60-
sigma.inv <- solve(Z %*% model$omega %*% t(Z) +
61-
diag(model$sigma2, nrow(Z)))
62-
6347
# calculate random effect by EB estimate
64-
re <- model$omega %*% t(Z) %*% sigma.inv %*% (y - X %*% beta)
48+
R <- solve(Z %*% model$omega %*% t(Z) + diag(model$sigma2, nrow(Z)),
49+
y - X %*% beta)
50+
re <- model$omega %*% t(Z) %*% R
6551

6652
# calculate broken stick estimate by summing fixed and random parts
6753
if (BS) re <- model$beta + re

R/brokenstick-class.R

Lines changed: 74 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -4,42 +4,93 @@
44
#' `brokenstick`. This object collects the fitted broken stick model.
55
#'
66
#' The package exports S3 methods for the `brokenstick` class for the following
7-
#' generic functions: [fitted()], [model.frame()], [model.matrix()], [plot()],
8-
#' [predict()], [print()], [residuals()].
7+
#' generic functions: [coef()], [fitted()], [model.frame()], [model.matrix()],
8+
#' [plot()], [predict()], [print()], [residuals()] and [summary()].
99
#'
10-
#' The package documents methods [fitted.brokenstick()], [plot.brokenstick()],
11-
#' [predict.brokenstick()] and [residuals.brokenstick()]. The package exports
12-
#' two helper functions for `brokenstick` objects: [get_knots()] and [get_r2()].
10+
#' The package exports the following helper functions for `brokenstick` objects:
11+
#' [get_knots()], [get_omega()] and [get_r2()].
1312
#'
1413
#' A `brokenstick` object is a list with the following named elements:
1514
#'
1615
#' @section Elements:
1716
#' \describe{
18-
#' \item{`names`}{A list with elements named `x`, `y` and `g` providing the
19-
#' variables names for the time, outcome and subject columns, respectively.}
20-
#' \item{`knots`}{Numeric vector of with the internal knots. Use [get_knots()]
21-
#' to extract knots.}
22-
#' \item{`boundary`}{Numeric vector of length 2 with the boundary knots. Use
23-
#' [get_knots()] to extract knots.}
24-
#' \item{`degree`}{The `degree` of the B-spline. See [splines::bs()]. Either
25-
#' 0 (constant model) or 1 (broken stick model).}
26-
#' \item{`model`}{Not used.}
27-
#' \item{`method`}{Either `lmer` or `kr`, identifying the fitting model.}
28-
#' \item{`control`}{List of control option returned by [control_brokenstick()].}
17+
#' \item{`call`}{Call that created the object}
18+
#' \item{`formula`}{A formula with the model specification,
19+
#' e.g. `formula(y ~ x | group)`}
20+
#' \item{`names`}{A named list with three elements (`"x"`, `"y"`, `"g"`)
21+
#' providing the variable name for time, outcome and subject, respectively.}
22+
#' \item{`internal`}{Numeric vector of with internal knots.}
23+
#' \item{`boundary`}{Numeric vector of length 2 with the boundary knots.}
24+
#' \item{`degree`}{The `degree` of the B-spline. See [splines::bs()]. Support
25+
#' only the values of 0 (step model) or 1 (broken stick model).}
26+
#' \item{`method`}{String, either `"kr"` or `"lmer"`, identifying the fitting model.}
27+
#' \item{`control`}{List of control options returned by [set_control()] used
28+
#' to set algorithmic details.}
2929
#' \item{`beta`}{Numeric vector with fixed effect estimates.}
3030
#' \item{`omega`}{Numeric matrix with variance-covariance estimates of the
31-
#' random effect.}
31+
#' broken stick estimates.}
3232
#' \item{`sigma2j`}{Numeric vector with estimates of the residual variance per
33-
#' group. Only used by method `kr`.}
33+
#' group. Only used by method `"kr"`.}
3434
#' \item{`sigma2`}{Numeric scalar with the mean residual variance.}
35-
#' \item{`draws`}{Numeric matrix with multiple imputations. The number of
36-
#' rows is equal to the number of missing values in `y`. The number of columns
37-
#' depends on `imp_skip`. Only used by `kr` if `imp_skip` is set.}
35+
#' \item{`sample`}{A numeric vector with descriptives of the training data.}
36+
#' \item{`light`}{Should the returned object be lighter? If `light = TRUE`
37+
#' the returned object will contain only the model settings and parameter
38+
#' estimates and not store the `sample`, `data`, `imp` and `mod` elements.
39+
#' The light object can be used to predict broken stick estimates for
40+
#' new data, but does not disclose the training data and is small.}
41+
#' \item{`data`}{The training data used to fit the model.}
42+
#' \item{`imp`}{The imputations generated for the missing outcome data. Only
43+
#' for `method = "kr"`.}
44+
#' \item{`mod`}{For `method = "kr"`: A named list with four components, each of
45+
#' class [coda::mcmc]. For `method = "lmer"`: An object of class
46+
#' [lme4::merMod-class].}
3847
#' }
3948
#'
4049
#' @name brokenstick-class
4150
#' @rdname brokenstick-class
42-
#' @author Stef van Buuren, 2020
43-
#' @references van Buuren S (2020).
51+
#' @author Stef van Buuren, 2021
52+
#' @references van Buuren S (2021).
4453
#' Broken Stick Model for Irregular Longitudinal Data. \emph{In preparation}.
4554
NULL
55+
56+
new_brokenstick <- function(call = match.call(),
57+
formula = formula(),
58+
names = list(x = character(),
59+
y = character(),
60+
g = character()),
61+
internal = numeric(0),
62+
boundary = numeric(0),
63+
degree = 1L,
64+
method = NA_character_,
65+
control = list(),
66+
beta = numeric(0),
67+
omega = numeric(0),
68+
sigma2j = numeric(0),
69+
sigma2 = numeric(0),
70+
light = FALSE,
71+
sample = numeric(0),
72+
data = numeric(0),
73+
imp = numeric(0),
74+
mod = list()) {
75+
result <- list(call = call,
76+
formula = formula,
77+
names = names,
78+
internal = internal,
79+
boundary = boundary,
80+
degree = degree,
81+
method = method,
82+
control = control,
83+
beta = beta,
84+
omega = omega,
85+
sigma2j = sigma2j,
86+
sigma2 = sigma2,
87+
sample = sample,
88+
light = light)
89+
if (!light) {
90+
result$data <- data
91+
result$imp <- imp
92+
result$mod <- mod
93+
}
94+
class(result) <- "brokenstick"
95+
return(result)
96+
}

R/brokenstick-constructor.R

Lines changed: 0 additions & 29 deletions
This file was deleted.

R/brokenstick-package.R

Lines changed: 9 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -13,33 +13,29 @@
1313
#' The main functions are:
1414
#' \tabular{ll}{
1515
#' \code{brokenstick()} \tab Fit a broken stick model to irregular data\cr
16-
#' \code{predict()} \tab Obtain predictions on new data\cr
1716
#' \code{plot()} \tab Plot observed and fitted trajectories by group \cr
17+
#' \code{predict()} \tab Obtain predictions on new data\cr
18+
#' \code{summary()} \tab Extract object summaries\cr
1819
#' }
1920
#'
2021
#' The following functions are user-oriented helpers:
2122
#' \tabular{ll}{
23+
#' \code{coef()} \tab Extract estimated parameters\cr
2224
#' \code{fitted()} \tab Calculate fitted values\cr
2325
#' \code{get_knots()} \tab Obtain the knots from a broken stick model\cr
26+
#' \code{get_omega()} \tab Extract variance-covariance of random effects\cr
2427
#' \code{get_r2()} \tab Obtain proportion of explained variance \cr
28+
#' \code{model.frame()} \tab Extract model frame\cr
29+
#' \code{model.matrix()} \tab Extract design matrix\cr
2530
#' \code{residuals()} \tab Extract residuals from broken stick model\cr
2631
#' }
2732
#'
28-
#' The following functions perform the calculations:
33+
#' The following functions perform calculations:
2934
#' \tabular{ll}{
30-
#' \code{control_brokenstick()}\tab Set controls to steer calculations\cr
31-
#' \code{EB()} \tab Empirical Bayes predictor for random effects\cr
32-
#' \code{kr()} \tab Kasim-Raudenbush sampler for two-level normal model \cr
33-
#' \code{make_basis()} \tab Create linear splines basis\cr
35+
#' \code{set_control()}\tab Set controls to steer calculations\cr
36+
#' \code{control_kr()}\tab Set controls for the \code{kr} method\cr
3437
#' }
3538
#'
36-
#' The package follows the \code{tidymodels} conventions
37-
#' \url{https://tidymodels.github.io/model-implementation-principles/}.
38-
#' For example, training data are not stored in the modelling object and
39-
#' calculated variables are named after the convention. The
40-
#' package architecture borrows important ideas from the \code{hardhat}
41-
#' package.(Vaughan, 2020)
42-
#'
4339
#' @docType package
4440
#' @name brokenstick-pkg
4541
#' @seealso \code{\link{brokenstick}},
@@ -51,7 +47,4 @@
5147
#' @references
5248
#' van Buuren, S. (2018). \emph{Flexible Imputation of Missing Data. Second Edition}. Chapman & Hall/CRC. Chapter 11.
5349
#' \url{https://stefvanbuuren.name/fimd/sec-rastering.html#sec:brokenstick}
54-
#'
55-
#' Vaughan, D. and Kuhn, M. (2020). \emph{hardhat: Construct Modeling Packages}.
56-
#' R package version 0.1.4. \url{https://CRAN.R-project.org/package=hardhat}
5750
NULL

0 commit comments

Comments
 (0)