Changelog

skmle 0.1.0

Bandwidth selection for the asynchronous estimators

kee_async(correct = "pairs") fits the pair-moment equation. Each subject’s own covariate pairs stand in for the covariate at the response time, reweighted by the ratio of response and covariate visit intensities estimated from the observation times. The equation has mean zero at the true coefficients for every bandwidth, with no model for the covariate process (no normality, no stationarity), and takes any number of covariates. correct = "gaussian" is the conditional-expectation correction for one stationary Gaussian covariate.
kee_async_cv() now selects the bandwidth for the corrected fit and returns the corrected fit. The default, method = "variance", selects the h minimising tr{Sigma_X V(h)}, with V(h) the sandwich variance of the corrected fit on all subjects and Sigma_X the second moment of X at the response times. The corrected fit is unbiased at every h, so this is its mean squared error; the selected h converges to the minimiser of the true variance and the Wald interval at it has its nominal level asymptotically. The default grid starts where at least ceiling(sqrt(n)) subjects have a pair inside the kernel.
method = "sync" selects the same fit’s bandwidth by synchronous-loss cross-validation. Held-out pair losses are regressed on the estimated variogram matrix at the pair’s lag and the intercept is minimised; a candidate whose intercept falls below the saturated loss on a fold is dropped. The default grid starts where every training fold has at least ceiling(sqrt(n_train)) subjects inside the kernel. The selected bandwidth gives a consistent estimate but does not converge, and the estimate can carry a selection bias of the order of its standard error. The previous criterion is available as method = "plain".
kee_async_cv() builds its default grid from the data. The candidates are n_h (default 100) log-spaced bandwidths from the smallest to the largest lag between a response and a covariate observation of the same subject (positive lags only when one_sided = TRUE). The previous default, log-spaced over 2(Q3 - Q1) n^-0.7 to 2(Q3 - Q1) n^-0.3 of the pooled observation times, had no justification for this estimator and is removed.
kee_async_cv() is plain K-fold cross-validation over subjects: each candidate is scored on the held-out subjects by the kernel-weighted squared error at that bandwidth, divided by the held-out weight, and the minimiser is selected. The new repeats argument averages the loss over several random splits (default 1).

Bias-corrected bandwidth selection

New functions skmle_bc(), kee_cox_bc() and kee_async_bc() choose the bandwidth and extrapolate the coefficients to h = 0. The three share one engine. Each fits its estimator (skmle(), kee_cox(), kee_async()) over a grid of bandwidths (by default 12, log-spaced) and chooses h by a plug-in estimate of the mean squared error: a weighted least-squares regression of the estimate on 1, h, [h^2,] 1/h across the grid, with weights from the sandwich variances. The coefficients are then extrapolated to h = 0 by a weighted regression on h over the grid points in [h/2, 2h]. The argument degree is 1 (the default) or 2; degree = 2 also removes an h^2 term. The returned objects have class "bc", with print() and summary() methods.
The reported covariance includes the choice of bandwidth. It is computed from the joint covariance of the estimates across bandwidths, built from each fit’s per-subject influence contributions, with a parametric bootstrap of the grid (n_boot, default 2000) that accounts for the data-driven choice of h. n_boot = 0 gives the variance with the bandwidth and window held fixed.
A Wald interval at the selected bandwidth undercovers. At the MSE-optimal bandwidth the bias and the standard deviation are of the same order. With bias of order h, the asymptotic coverage of the nominal 95% interval is 0.891.
In simulation, the 95% intervals for the first coefficient from the new functions had coverage between 0.929 and 0.960. With 1000 replicates per cell and n from 200 to 3200: for skmle() with s in {0, 0.5, 1}, the coverage of the interval for beta_1 was 0.936 to 0.960 with skmle_bc(), against 0.720 to 0.770 at n = 3200 for the interval at the selected bandwidth without extrapolation. kee_cox_bc() gave 0.938 to 0.956. For the slope, kee_async_bc(), with an Ornstein-Uhlenbeck covariate (covariance exp(-|s - t|)) and a two-sided kernel, gave 0.929 to 0.950, while the split-half selector of Cao, Zeng and Fine (2015) gave 0.869 at n = 200, falling to 0.783 at n = 3200.
skmle(), kee_cox() and kee_async() return scores. These are the per-subject contributions to the score; their crossproduct reproduces the meat of the sandwich.
?skmle_cv is corrected: the carried-forward approximation in the held-out score does not cancel across bandwidths. The help page said it did. The held-out score rewards fits attenuated by the time since the last visit, so the selected bandwidth converges to a constant fixed by the visit process rather than to zero, and the coefficients at h_cv are biased for beta. In simulation (200 replicates per cell, n = 100 to 400, grid [0.03, 1]), the coverage of the interval for beta_1 at the cross-validated bandwidth was 0.760 to 0.880. Use skmle_cv() for prediction and skmle_bc() for inference on beta.

Bug fixes

kee_async(correct = "pairs") no longer fails when the response and covariate schedules separate in time. The response/covariate intensity ratio comes from a logistic fit on the observation times; with a covariate collected in one late window plus a stray early draw, the fitted probability reached 0 or 1, the ratio became Inf or NaN, and the pair fit stopped with “missing value where TRUE/FALSE needed”. The probability is now bounded at 0.5/N from either end, with a warning when the bound binds, and an aliased spline coefficient is taken as zero.
skmle_cv() scores a failed or overflowing fold as +Inf, and refuses a grid with no finite candidate. A held-out score whose linear predictor overflowed returned Inf - Inf = NaN, which propagated into the candidate’s mean loss; which.min() skipped it, and a grid on which every candidate was NaN crashed on an NA condition. Worse, a grid on which every candidate was +Inf (every training fit failed) selected its first entry. Both cases now count as a lost candidate, a grid with none finite stops with a message, and quiet = FALSE reports how many candidates were ineligible. Finite losses are unchanged to the bit. Found on ADNI MCI-to-dementia data, where the smaller default-grid bandwidths overflowed on some folds.
A custom weight now reaches the risk-set averages. kee_cox(), kee_additive() and the skmle() sandwich used the weight only for the row weights and still built every risk-set average S^(k) from a hardcoded Epanechnikov kernel. A weight other than the default therefore gave a kee_cox() or kee_additive() estimate solving an equation that belongs to no estimator, and skmle() standard errors for a different estimator from the one fitted. The default kernels are unaffected: their estimates are bit-identical to before.
one_sided = FALSE now reaches the whole variance. calc_B() (the skmle() meat) and the risk set in kee_additive()’s B required a positive lag whatever one_sided said. The full-kernel scheme drops the indicator in every term (Sun et al. 2025, Remark 7; Cao et al. 2015, eqs. 3.8-3.9). With the full kernel, skmle() standard errors and kee_additive() estimates change.
The skmle() sandwich evaluates the baseline at the risk-set time. The risk-set averages in calc_A() and calc_B() used alpha(X_k), the baseline at each risk-set member’s own event time, where Theorem 3.5 and the R prototype use alpha(X_i), the time the risk set is formed at. Standard errors change for every s except 0.5, where H''/H = 1 and the two agree. This had been present since the first commit.
The skmle() meat sums each subject’s cumulative-hazard term. The second term of the meat averaged over a subject’s observations, following the R prototype, where the published Omega-hat and the objective’s own per-subject score sum over them. In a simulation with 2,000 replicates per cell (s = 0, 0.5, 1; n = 200, 400), the averaged form gave SE/ESE of 0.90 to 0.96 and coverage of 91.5% to 94.2%; the summed form gave 0.93 to 0.98 and 92.5% to 95.0%, closer to 95% in all 12 comparisons.
skmle() no longer returns a stalled optimisation. With s != 0, SLSQP occasionally stalled on the floor of the cumulative-hazard integrand (0.3% of fits at s = 1 in simulation, nlopt status -1), or stopped with an event row on the feasibility boundary under a success status (0.4%), and skmle() returned those fits silently with standard errors between 1e-16 and 1e-3. Such a fit is now restarted: the problem is solved once with the quadrature pairs also constrained, and the published problem is re-solved from that point, in skmle() and in the cross-validation loop alike. The restart is kept only if it reaches at least as good an objective. The estimator is unchanged; constraining the quadrature pairs outright removes the stalls but biased beta by up to 0.35 empirical SDs at s = 1, n = 200. skmle() warns if the optimiser still reports failure.
Rows with a negative weight are kept by the cumulative-hazard term of the objective and by the meat, as they already were by the event term and the bread.
For s != 0 the feasibility constraints cover every row the weight reaches, not only lags within h. For the Epanechnikov kernels the set is unchanged.
weight now comes after one_sided in skmle(), kee_cox(), kee_additive() and skmle_cv(), so positional calls written before weight existed work again.

Any weight function, not only the ones supplied

skmle(), kee_cox(), kee_additive() and skmle_cv() take a weight argument. Any R function of the standardised lag u = (t - r)/h may be passed; the fitting code privileges no particular kernel. w_epan_half() and w_epan_full() build the two Epanechnikov kernels and are exported.
A weight takes one argument and returns the shape it was given. The sieve quadrature evaluates it on a matrix of node-by-observation lags, not only on a vector, so a weight built from apply() or a scalar if() would be recycled into the wrong cells. Both the length and the matrix shape are checked when the argument is resolved, rather than failing later inside the fit.
A weight may take negative values and need not integrate to 1. The normalisation cancels in the risk-set ratios and is absorbed into h elsewhere, which is why w_epan_half() integrates to 1/2.
weight defaults to NULL, meaning the Epanechnikov that one_sided implies. A literal w_epan_half() default would have been wrong, and was wrong for one commit: its support is [0, 1], so one_sided = FALSE with that default zeroes every negative lag through the support and returns a one-sided fit under a two-sided name. The package’s own suite caught it and now pins it.
one_sided stays a separate argument, because it is a property of the model and not of the kernel. It is the risk-set restriction of a hazard model and applies to the sign of the lag after the weight is evaluated, so a two-sided weight under one_sided = TRUE is truncated at zero rather than rejected.
kee_async(), kee_async_td() and kee_async_cv() do not take weight; they use the full Epanechnikov of Cao, Zeng and Fine.

The cross-validation loop runs in C++

The fold loop and the held-out score moved back into C++, as skmle_cv_cpp(). The reason is not speed, which is dominated by the fit either way. Scoring in R required a second implementation of the model’s transformation: R/utils.R carried a trans_link() that mirrored trans_fun() in src/skmle_cpp.cpp by hand, floor and all, and because trans_fun is not exported nothing could compare them and no test did. The criterion and the objective now call the same function, which is a guarantee rather than a test someone has to remember to write. trans_link() is gone.
The optimiser set-up is factored into one internal fit_core(), shared by skmle_cpp_fit() and the loop, so the fit the cross-validation scores cannot drift from the fit skmle() returns.
locf_score() stays in R. It is set-up rather than a loop, it needs splines::ns() and findInterval(), and it evaluates no part of the model: it returns a basis, quadrature weights and row indices.
The criterion is unchanged, and that is checked rather than asserted. The losses from the R loop were recorded before the port, and the C++ loop reproduces them to 4.4e-16 with identical coefficients. test-cv-reference.R keeps the comparison in the suite.
The fold loop evaluates the caller’s weight function itself, as an Rcpp::Function, so the kernel keeps one description. An earlier draft of this change gave C++ a polynomial description of the weight alongside the R closure and a test comparing the two; handing the closure across removes the second description rather than testing it.

Covariate observation times may be negative

obs_times is no longer required to be non-negative. Event times still are, because the sieve basis and the cumulative-hazard quadrature are built on [0, max(X)] and anything outside that span is extrapolation. An observation time is never an argument to a basis: it reaches the score only through the lag X - obs_times, and the one-sided rule tests the sign of that lag, not the sign of the time. The old check applied a constraint belonging to one argument to an argument that does not carry it.
What this buys is a covariate measured before the participant entered follow-up, a screening draw or a run-in visit, without re-origining the time axis. Shifting the origin moves X too and leaves every lag and every risk set unchanged, so it was never a way to recover that information.
The permission is numerical and the interpretation is not. A pre-entry measurement has to be commensurable with the on-study ones, entry may not be independent of the covariate that preceded it, and a participant contributes a pre-entry row only by having survived to enrol. None of this is checked. ?kee_cox states it under “Negative observation times”.
Reach is unchanged: a row further back than one bandwidth from every event time still contributes nothing, and the fit equals the fit with that row deleted.

Weights that can take negative values

The estimating equations and the sieve likelihood guard on a nonzero row weight rather than a positive one. A weight function supplied by a caller is not required to be nonnegative, and the old > 0 test dropped negative-weight rows silently instead of failing, which made the score wrong with no diagnostic. The objective, its gradient, and the bread and meat of the sandwich now sum over the same rows; previously a negative weight would have put them on different row sets.
The internal risk-set weights keep their > 0 test, which is a skip-zero optimisation on a kernel that is nonnegative by construction, not a sign assumption about the caller’s weight.
The kernels shipped with the package are nonnegative, so every result the package has produced is unchanged to the last bit; the difference appears only under a caller-supplied weight.

Scope

The package covers two settings, not one. Alongside the transformed hazards models for survival outcomes it now fits generalised linear models for asynchronous longitudinal outcomes, where the response and the covariate are recorded on different time grids. The Title drops its “for Survival Models” restriction, and the Description, README, package help page and tutorial have been rewritten to present the two settings on equal footing rather than treating the second as an extra.

Getting a first answer

The package assumed you already knew the method. Three things a newcomer could not supply now have defaults, each announced rather than hidden.

h is optional in every fitting function. When omitted, a rule-of-thumb bandwidth is read off the observation times – the geometric midpoint of the grid the matching _cv() function searches – and a message reports the value, says it is a starting point rather than a tuned choice, and names the function that tunes it. Supplying h keeps the message quiet.
s defaults to 0 in skmle() and skmle_cv(), the proportional hazards model, so the Box-Cox parameter no longer has to be understood before the first fit.
times defaults in kee_async_td() to 25 points spanning the 10th to 90th percentile of the observed response times, which keeps them away from the edges where a one-sided window makes the curve unreliable.

Together these mean kee_cox(Surv(X, delta) ~ z, data, id = id, obs_times = obs_times) and data_y |> kee_async(data_x, y ~ x, id = id, time = time) are complete calls.

print() and summary() now name the model in words and report the bandwidth and kernel, so the output says what was fitted rather than assuming you remember.
plot() on a cv.skmle or cv.kee_async object draws the held-out loss against bandwidth with the selection marked, which is the quickest way to see a minimum sitting on the edge of the grid.
A “Which function do I need?” table in ?skmle-package, the README and the tutorial.

Interface

The asynchronous estimators take the two data frames first, so they compose with the native pipe: data_y |> kee_async(data_x, y ~ x, id = id, time = time, h = 0.25). The formula still spans both tables – its left-hand side is looked up in data_y, its right-hand side in data_x – which is unavoidable when the data genuinely lives in two frames, and is now stated in the argument docs rather than left to be discovered.
id and time accept a bare name, a string, or {{ col }}, so the estimators can be wrapped in other functions. They previously used substitute(), which deparsed {{ col }} literally.
tidy(), glance() and augment() methods. augment() attaches .fitted to the covariate table, which is where a fitted value lives here; there is no .resid, because a residual would need a response value at the covariate time and that is exactly what asynchronous data lacks.
sim_async_data(), confint(), and both cv_results tables return tibbles.
Requires R (>= 4.1) for the native pipe used throughout the examples.

Usability

kee_async() and kee_async_td() take a formula spanning the two tables – kee_async(y ~ x, data_y, data_x, id, time, h) – instead of seven positional vectors. Swapping the response and covariate tables used to return plausible numbers with no complaint.
kee_async_cv() selects the bandwidth by subject-level cross-validation, scoring candidates by kernel-weighted squared error on held-out subjects. The asynchronous estimators previously offered no guidance on h at all.
kee_cox() and kee_additive() now check that times lie on [0, 1], which skmle() already did. kee_additive() genuinely requires it – its quadrature is built on [0, 1] – so times in other units were silently wrong rather than merely unusual.
kee_async() warns when fewer than 5% of response occasions have a covariate observation in their window, which is the signature of a bandwidth on the wrong scale. The asynchronous estimators are scale-free, so this cannot be checked by a range test.
vcov(), nobs() and (for kee_td) confint() methods. confint() previously failed on every fitted object in the package, because there was no vcov() for the default method to call.
skmle() and skmle_cv() no longer take norder. It was validated and documented but never used: the sieve basis is a natural cubic spline, whose order is fixed. Existing calls that pass it will now error, and should drop the argument.
sim_async_data() returns the covariates as plain columns (x, or x1, x2, …) rather than a matrix column, so they can be named in a formula.
New article, vignette("asynchronous"): why last-value-carried-forward and regression calibration are inconsistent here, how to read the bandwidth sensitivity plot, what the half kernel changes, and how to get the units right.

Asynchronous longitudinal data

kee_async() and kee_async_td() implement the kernel-weighted estimating equations of Cao, Zeng and Fine (2015) for a longitudinal response and a longitudinal covariate observed on different time grids, with time-invariant and time-dependent coefficients respectively. Identity, log and logistic links are supported.
sim_async_data() simulates from their Section 4 design.
plot() and print() methods for the kee_td coefficient curves.
Both estimators are backed by C++: pairs are enumerated once and collapsed onto covariate rows, so the pair design never enters the Newton loop. The time-dependent weight factorises over the two occasion indices, so no pair is enumerated there at all.

Bandwidth selection in `skmle_cv()`

The held-out loss no longer falls away with the bandwidth. It was the kernel-weighted log-likelihood evaluated on the held-out fold, and the kernel weight is W(u/h)/h, so the weight a subject contributes shrinks as h grows and the criterion decreased monotonically whatever the fit was worth. The largest candidate in h_grid won every grid on every data set – including the automatic grid – and coefficients further from the truth were preferred to coefficients nearer it. Reported by a reader of the tutorial, whose CV curve was monotone across n = 200/1000, censoring 0.2/0.7 and five seeds.
The score is now an ordinary log-likelihood evaluated on the held-out subjects, with the covariate path carried forward from the last observation. No kernel and no bandwidth enter it, so it is the same yardstick for every candidate and the comparison is about the fit. See ?skmle_cv.
A minimum at an endpoint of h_grid now warns, as kee_async_cv() already did: the selected value is then the best of the values offered rather than a minimum. On the automatic grid this fires often, because the grid stops at the rate-based tau * n^-0.3 while the finite-sample minimum of the loss commonly lies above it.
cv_results gains an se column, the standard error of each loss across the folds, because the criterion is flat: at the sample sizes in the tutorial every candidate is within one standard error of the minimum. It also scores prediction of the held-out hazard, where the baseline can absorb attenuation in beta, so it leans towards more smoothing than the coefficients on their own would want. Over 10 replicates at n = 200 on a grid spanning 0.05 to 0.9, the squared error of beta-hat at the selected bandwidth averaged 0.155, against 0.205 at the largest candidate – what the old criterion always returned – and 0.065 at the oracle bandwidth. A one-standard-error rule does not help: it lands in the noisy small-h end and scored 0.205, and ?skmle_cv says so.
h_grid = NULL worked again. The automatic grid referred to the end of follow-up before it was computed, so the documented default errored with object 'tau' not found for any data set.
The fold loop moved from C++ to R. It only subsets matrices and calls the optimiser, which R does in a dozen lines, and the held-out score needs a spline basis at subject-specific quadrature nodes.

Half and full kernels

skmle(), kee_cox(), kee_additive(), skmle_cv(), kee_async() and kee_async_td() take a one_sided argument. TRUE is the half kernel: only covariate observations preceding a time inform it. FALSE is the full, two-sided kernel. The survival estimators default to TRUE, the estimator as published; the asynchronous ones default to FALSE, the kernel of Cao, Zeng and Fine.
For the survival estimators the switch reaches the risk-set averages in the C++ backend as well as the row weights, so the two halves of an estimator cannot disagree. skmle_cv() also uses it inside the fold loop, which previously had the half kernel hardcoded, so the bandwidth is now selected under the kernel the refit uses.
The Epanechnikov kernel had been written out inline in three fitting functions; it now lives in one internal helper, so the half/full switch cannot be applied to two of the three.

Initial release

Initial CRAN release.
skmle() fits transformed hazards models by sieve maximum kernel-weighted log-likelihood estimation (SMKLE).
kee_cox() and kee_additive() fit Cox and additive hazards models via kernel-weighted estimating equations.
skmle_cv() selects the kernel bandwidth by K-fold cross-validation.
sim_skmle_data() simulates survival data with sparse, intermittently observed longitudinal covariates.
plot(), print() and summary() methods for fitted objects.