For right-censored data, fit a regularized Cox cure rate model through elastic-net penalty following Masud et al. (2018), and Zou and Hastie (2005). For right-censored data with uncertain event status, fit the regularized Cox cure model proposed by Wang et al. (2020). Without regularization, the model reduces to the regular Cox cure rate model (Kuk and Chen, 1992; Sy and Taylor, 2000)

cox_cure_net(
  surv_formula,
  cure_formula,
  time,
  event,
  data,
  subset,
  contrasts = NULL,
  surv_lambda = NULL,
  surv_alpha = 1,
  surv_nlambda = 10,
  surv_lambda_min_ratio = 0.1,
  surv_l1_penalty_factor = NULL,
  cure_lambda = NULL,
  cure_alpha = 1,
  cure_nlambda = 10,
  cure_lambda_min_ratio = 0.1,
  cure_l1_penalty_factor = NULL,
  cv_nfolds = 0,
  surv_start = NULL,
  cure_start = NULL,
  surv_offset = NULL,
  cure_offset = NULL,
  surv_standardize = TRUE,
  cure_standardize = TRUE,
  em_max_iter = 200,
  em_rel_tol = 1e-05,
  surv_max_iter = 10,
  surv_rel_tol = 1e-05,
  cure_max_iter = 10,
  cure_rel_tol = 1e-05,
  tail_completion = c("zero", "exp", "zero-tau"),
  tail_tau = NULL,
  pmin = 1e-05,
  early_stop = TRUE,
  verbose = FALSE,
  ...
)

cox_cure_net.fit(
  surv_x,
  cure_x,
  time,
  event,
  cure_intercept = TRUE,
  surv_lambda = NULL,
  surv_alpha = 1,
  surv_nlambda = 10,
  surv_lambda_min_ratio = 0.1,
  surv_l1_penalty_factor = NULL,
  cure_lambda = NULL,
  cure_alpha = 1,
  cure_nlambda = 10,
  cure_lambda_min_ratio = 0.1,
  cure_l1_penalty_factor = NULL,
  cv_nfolds = 0,
  surv_start = NULL,
  cure_start = NULL,
  surv_offset = NULL,
  cure_offset = NULL,
  surv_standardize = TRUE,
  cure_standardize = TRUE,
  em_max_iter = 200,
  em_rel_tol = 1e-05,
  surv_max_iter = 10,
  surv_rel_tol = 1e-05,
  cure_max_iter = 10,
  cure_rel_tol = 1e-05,
  tail_completion = c("zero", "exp", "zero-tau"),
  tail_tau = NULL,
  pmin = 1e-05,
  early_stop = TRUE,
  verbose = FALSE,
  ...
)

Arguments

surv_formula

A formula object starting with ~ for the model formula in survival model part. For Cox model, no intercept term is included even if an intercept is specified or implied in the model formula. A model formula with an intercept term only is not allowed.

cure_formula

A formula object starting with ~ for the model formula in incidence model part. For logistic model, an intercept term is included by default and can be excluded by adding + 0 or - 1 to the model formula.

time

A numeric vector for the observed survival times.

event

A numeric vector for the event indicators. NA's are allowed and represent uncertain event status.

data

An optional data frame, list, or environment that contains the covariates and response variables (time and event), included in the model. If they are not found in data, the variables are taken from the environment of the specified formula, usually the environment from which this function is called.

subset

An optional logical vector specifying a subset of observations to be used in the fitting process.

contrasts

An optional list, whose entries are values (numeric matrices or character strings naming functions) to be used as replacement values for the contrasts replacement function and whose names are the names of columns of data containing factors. See contrasts.arg of model.matrix.default for details.

surv_lambda, cure_lambda

A numeric vector consists of nonnegative values representing the tuning parameter sequence for the survival model part or the incidence model part.

surv_alpha, cure_alpha

A number between 0 and 1 for tuning the elastic net penalty for the survival model part or the incidence model part. If it is one, the elastic penalty will reduce to the well-known lasso penalty. If it is zero, the ridge penalty will be used.

surv_nlambda, cure_nlambda

A positive number specifying the number of surv_lambda or cure_lambda if surv_lambda or cure_lambda is not specified, respectively. The default value is 10.

surv_lambda_min_ratio, cure_lambda_min_ratio

The ratio of the minimum surv_lambda (or cure_lambda) to the large enough surv_lambda (or codecure_lambda) that produces all-zero estimates on log scale. The default value is 1e-1.

surv_l1_penalty_factor, cure_l1_penalty_factor

A numeric vector that consists of nonnegative penalty factors (or weights) on L1-norm for the coefficient estimate vector in the survival model part or the incidence model part. The penalty is applied to the coefficient estimate divided by the specified weights. The specified weights are re-scaled internally so that their summation equals the length of coefficients. If NULL is specified, the weights are all set to be one.

cv_nfolds

An non-negative integer specifying number of folds in cross-validation (CV). The default value is 0 and the CV procedure is not enabled.

surv_start

An optional numeric vector representing the starting values for the survival model component or the incidence model component. If surv_start = NULL is specified, the starting values will be obtained from fitting a regular Cox to events only. Similarly, if cure_start = NULL is specified, the starting values will be obtained from fitting a regular logistic model to the non-missing event indicators.

cure_start

An optional numeric vector representing the starting values for the survival model component or the incidence model component. If surv_start = NULL is specified, the starting values will be obtained from fitting a regular Cox to events only. Similarly, if cure_start = NULL is specified, the starting values will be obtained from fitting a regular logistic model to the non-missing event indicators.

surv_offset

An optional numeric vector representing the offset term in the survival model compoent or the incidence model component. The function will internally try to find values of the specified variable in the data first. Alternatively, one or more offset terms can be specified in the formula (by stats::offset()). If more than one offset terms are specified, their sum will be used.

cure_offset

An optional numeric vector representing the offset term in the survival model compoent or the incidence model component. The function will internally try to find values of the specified variable in the data first. Alternatively, one or more offset terms can be specified in the formula (by stats::offset()). If more than one offset terms are specified, their sum will be used.

surv_standardize, cure_standardize

A logical value specifying whether to standardize the covariates for the survival model part or the incidence model part. If FALSE, the covariates will be standardized internally to have mean zero and standard deviation one.

em_max_iter

A positive integer specifying the maximum iteration number of the EM algorithm. The default value is 200.

em_rel_tol

A positive number specifying the tolerance that determines the convergence of the EM algorithm in terms of the convergence of the covariate coefficient estimates. The tolerance is compared with the relative change between estimates from two consecutive iterations, which is measured by ratio of the L1-norm of their difference to the sum of their L1-norm. The default value is 1e-5.

surv_max_iter, cure_max_iter

A positive integer specifying the maximum iteration number of the M-step routine related to the survival model component or the incidence model component. The default value is 10 to encourage faster convergence.

surv_rel_tol

A positive number specifying the tolerance that determines the convergence of the M-step related to the survival model component or the incidence model component in terms of the convergence of the covariate coefficient estimates. The tolerance is compared with the relative change between estimates from two consecutive iterations, which is measured by ratio of the L1-norm of their difference to the sum of their L1-norm. The default value is 1e-5.

cure_rel_tol

A positive number specifying the tolerance that determines the convergence of the M-step related to the survival model component or the incidence model component in terms of the convergence of the covariate coefficient estimates. The tolerance is compared with the relative change between estimates from two consecutive iterations, which is measured by ratio of the L1-norm of their difference to the sum of their L1-norm. The default value is 1e-5.

tail_completion

A character string specifying the tail completion method for conditional survival function. The available methods are "zero" for zero-tail completion after the largest event times (Sy and Taylor, 2000), "exp" for exponential-tail completion (Peng, 2003), and "zero-tau" for zero-tail completion after a specified tail_tau. The default method is the zero-tail completion proposed by Sy and Taylor (2000).

tail_tau

A numeric number specifying the time of zero-tail completion. It will be used only if tail_completion = "zero-tau". A reasonable choice must be a time point between the largest event time and the largest survival time.

pmin

A numeric number specifying the minimum value of probabilities for sake of numerical stability. The default value is 1e-5.

early_stop

A logical value specifying whether to stop the iteration once the negative log-likelihood unexpectedly increases, which may suggest convergence on likelihood, or indicate numerical issues or implementation bugs. The default value is TRUE.

verbose

A logical value. If TRUE, a verbose information will be given along iterations for tracing the convergence. The default value is FALSE.

...

Other arguments for future usage. A warning will be thrown if any invalid argument is specified.

surv_x

A numeric matrix for the design matrix of the survival model component.

cure_x

A numeric matrix for the design matrix of the cure rate model component. The design matrix should exclude an intercept term unless we want to fit a model only including the intercept term. In that case, we need further set cure_intercept = FALSE to not standardize the intercept term.

cure_intercept

A logical value specifying whether to add an intercept term to the cure rate model component. If TRUE by default, an intercept term is included.

Value

cox_cure_net object for regular Cox cure rate model or cox_cure_net_uncer object for Cox cure rate model with uncertain events.

Details

The model estimation procedure follows expectation maximization (EM) algorithm. Variable selection procedure through regularization by elastic net penalty is developed based on cyclic coordinate descent and majorization-minimization (MM) algorithm.

References

Kuk, A. Y. C., & Chen, C. (1992). A mixture model combining logistic regression with proportional hazards regression. Biometrika, 79(3), 531--541.

Masud, A., Tu, W., & Yu, Z. (2018). Variable selection for mixture and promotion time cure rate models. Statistical methods in medical research, 27(7), 2185--2199.

Peng, Y. (2003). Estimating baseline distribution in proportional hazards cure models. Computational Statistics & Data Analysis, 42(1-2), 187--201.

Sy, J. P., & Taylor, J. M. (2000). Estimation in a Cox proportional hazards cure model. Biometrics, 56(1), 227--236.

Wang, W., Luo, C., Aseltine, R. H., Wang, F., Yan, J., & Chen, K. (2020). Suicide Risk Modeling with Uncertain Diagnostic Records. arXiv preprint arXiv:2009.02597.

Zou, H., & Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2), 301--320.

See also

cox_cure for regular Cox cure rate model.

Examples

library(intsurv) ### regularized Cox cure rate model ================================== ## simulate a toy right-censored data with a cure fraction set.seed(123) n_obs <- 100 p <- 10 x_mat <- matrix(rnorm(n_obs * p), nrow = n_obs, ncol = p) colnames(x_mat) <- paste0("x", seq_len(p)) surv_beta <- c(rep(0, p - 5), rep(1, 5)) cure_beta <- c(rep(1, 2), rep(0, p - 2)) dat <- simData4cure(nSubject = n_obs, lambda_censor = 0.01, max_censor = 10, survMat = x_mat, survCoef = surv_beta, cureCoef = cure_beta, b0 = 0.5, p1 = 1, p2 = 1, p3 = 1) ## model-fitting from given design matrices fit1 <- cox_cure_net.fit(x_mat, x_mat, dat$obs_time, dat$obs_event, surv_nlambda = 10, cure_nlambda = 10, surv_alpha = 0.8, cure_alpha = 0.8) ## model-fitting from given model formula fm <- paste(paste0("x", seq_len(p)), collapse = " + ") surv_fm <- as.formula(sprintf("~ %s", fm)) cure_fm <- surv_fm fit2 <- cox_cure_net(surv_fm, cure_fm, data = dat, time = obs_time, event = obs_event, surv_alpha = 0.5, cure_alpha = 0.5) ## summary of BIC's BIC(fit1)
#> df surv_df cure_df BIC #> 1 2 0 2 571.6569 #> 2 2 0 2 567.3181 #> 3 3 0 3 566.3957 #> 4 4 0 4 567.0484 #> 5 6 0 6 572.1122 #> 6 7 0 7 572.3148 #> 7 9 0 9 578.0111 #> 8 10 0 10 579.5357 #> 9 10 0 10 577.3777 #> 10 11 0 11 580.3645 #> 11 3 1 2 574.8266 #> 12 3 1 2 570.4866 #> 13 4 1 3 569.5839 #> 14 5 1 4 570.2510 #> 15 6 1 5 570.7119 #> 16 8 1 7 575.5220 #> 17 10 1 9 581.2231 #> 18 11 1 10 582.7494 #> 19 11 1 10 580.5933 #> 20 12 1 11 583.5942 #> 21 4 2 2 571.2981 #> 22 4 2 2 566.9017 #> 23 5 2 3 565.9767 #> 24 6 2 4 566.6088 #> 25 7 2 5 567.0865 #> 26 9 2 7 571.9303 #> 27 11 2 9 577.6616 #> 28 12 2 10 579.1972 #> 29 12 2 10 577.0534 #> 30 13 2 11 580.0550 #> 31 6 4 2 572.0628 #> 32 6 4 2 567.2647 #> 33 7 4 3 566.1255 #> 34 8 4 4 566.5304 #> 35 9 4 5 567.1024 #> 36 11 4 7 572.0521 #> 37 12 4 8 573.2681 #> 38 14 4 10 579.5170 #> 39 14 4 10 577.4542 #> 40 15 4 11 580.4828 #> 41 7 5 2 561.6637 #> 42 7 5 2 556.4483 #> 43 8 5 3 555.4625 #> 44 9 5 4 555.9295 #> 45 9 5 4 552.1100 #> 46 12 5 7 561.6186 #> 47 13 5 8 562.8272 #> 48 15 5 10 569.1044 #> 49 15 5 10 567.0317 #> 50 16 5 11 570.0577 #> 51 8 6 2 553.9714 #> 52 8 6 2 549.1187 #> 53 9 6 3 548.2548 #> 54 10 6 4 548.7672 #> 55 10 6 4 544.9904 #> 56 13 6 7 554.6001 #> 57 15 6 9 560.3878 #> 58 16 6 10 562.0681 #> 59 16 6 10 559.9757 #> 60 17 6 11 562.9819 #> 61 8 6 2 546.0186 #> 62 9 7 2 545.7555 #> 63 10 7 3 544.8851 #> 64 11 7 4 545.4099 #> 65 11 7 4 541.6443 #> 66 14 7 7 551.3362 #> 67 16 7 9 557.1108 #> 68 17 7 10 558.8130 #> 69 17 7 10 556.7136 #> 70 18 7 11 559.7027 #> 71 10 8 2 549.1734 #> 72 10 8 2 544.4072 #> 73 11 8 3 543.5357 #> 74 12 8 4 544.0616 #> 75 12 8 4 540.3017 #> 76 14 8 6 545.4479 #> 77 17 8 9 555.8321 #> 78 18 8 10 557.5624 #> 79 18 8 10 555.4684 #> 80 19 8 11 558.4443 #> 81 12 10 2 553.4890 #> 82 12 10 2 548.5615 #> 83 13 10 3 547.6633 #> 84 14 10 4 548.1603 #> 85 14 10 4 544.3829 #> 86 16 10 6 549.5268 #> 87 19 10 9 559.9157 #> 88 20 10 10 561.6565 #> 89 20 10 10 559.5554 #> 90 21 10 11 562.5128 #> 91 12 10 2 549.7380 #> 92 12 10 2 545.0480 #> 93 13 10 3 544.0637 #> 94 14 10 4 544.5302 #> 95 14 10 4 540.7353 #> 96 16 10 6 545.8863 #> 97 19 10 9 556.3017 #> 98 20 10 10 558.0478 #> 99 21 10 11 560.5439 #> 100 21 10 11 558.8899
BIC(fit2)
#> df surv_df cure_df BIC #> 1 2 0 2 571.7958 #> 2 2 0 2 568.4844 #> 3 3 0 3 568.4511 #> 4 4 0 4 569.2552 #> 5 6 0 6 574.2389 #> 6 7 0 7 574.6851 #> 7 8 0 8 575.6070 #> 8 10 0 10 581.4934 #> 9 10 0 10 578.9973 #> 10 11 0 11 581.7157 #> 11 3 1 2 575.2438 #> 12 3 1 2 571.9320 #> 13 4 1 3 571.9106 #> 14 5 1 4 572.7245 #> 15 7 1 6 577.7091 #> 16 8 1 7 578.1582 #> 17 9 1 8 579.0867 #> 18 11 1 10 584.9727 #> 19 11 1 10 582.4786 #> 20 12 1 11 585.2040 #> 21 4 2 2 573.4404 #> 22 4 2 2 570.0984 #> 23 5 2 3 570.0641 #> 24 6 2 4 570.8632 #> 25 8 2 6 575.8681 #> 26 9 2 7 576.3247 #> 27 11 2 9 581.8671 #> 28 12 2 10 583.1623 #> 29 12 2 10 580.6766 #> 30 13 2 11 583.4063 #> 31 6 4 2 576.5390 #> 32 6 4 2 573.0304 #> 33 7 4 3 572.8699 #> 34 8 4 4 573.5255 #> 35 10 4 6 578.5651 #> 36 11 4 7 579.1056 #> 37 13 4 9 584.6872 #> 38 14 4 10 586.0446 #> 39 14 4 10 583.6066 #> 40 15 4 11 586.3663 #> 41 7 5 2 569.3109 #> 42 7 5 2 565.4851 #> 43 8 5 3 565.3667 #> 44 9 5 4 566.0206 #> 45 11 5 6 571.2470 #> 46 12 5 7 571.7412 #> 47 14 5 9 577.3128 #> 48 15 5 10 578.7157 #> 49 15 5 10 576.2806 #> 50 16 5 11 579.0345 #> 51 8 6 2 562.1769 #> 52 8 6 2 558.4644 #> 53 9 6 3 558.4153 #> 54 10 6 4 559.0812 #> 55 11 6 5 559.7931 #> 56 13 6 7 564.9129 #> 57 15 6 9 570.4631 #> 58 16 6 10 571.8919 #> 59 16 6 10 569.4491 #> 60 17 6 11 572.1968 #> 61 8 6 2 553.2210 #> 62 8 6 2 549.5238 #> 63 9 6 3 549.5231 #> 64 10 6 4 550.2024 #> 65 11 6 5 550.9437 #> 66 13 6 7 556.1250 #> 67 15 6 9 561.6600 #> 68 16 6 10 563.1094 #> 69 16 6 10 560.6603 #> 70 17 6 11 563.3988 #> 71 10 8 2 555.7105 #> 72 10 8 2 551.8308 #> 73 11 8 3 551.8037 #> 74 12 8 4 552.4865 #> 75 13 8 5 553.2438 #> 76 15 8 7 558.4816 #> 77 17 8 9 564.0158 #> 78 18 8 10 565.4922 #> 79 18 8 10 563.0476 #> 80 19 8 11 565.7790 #> 81 12 10 2 558.9853 #> 82 12 10 2 555.1329 #> 83 13 10 3 555.0467 #> 84 14 10 4 555.7139 #> 85 15 10 5 556.4813 #> 86 17 10 7 561.7638 #> 87 19 10 9 567.2879 #> 88 20 10 10 568.7817 #> 89 20 10 10 566.3330 #> 90 21 10 11 569.0486 #> 91 12 10 2 554.2454 #> 92 12 10 2 550.4021 #> 93 13 10 3 550.3402 #> 94 14 10 4 551.0016 #> 95 15 10 5 551.7540 #> 96 17 10 7 557.0672 #> 97 19 10 9 562.5825 #> 98 20 10 10 564.0952 #> 99 20 10 10 561.6475 #> 100 21 10 11 564.3541
## list of coefficient estimates based on BIC coef(fit1)
#> $surv #> x1 x2 x3 x4 x5 x6 #> 0.00000000 0.04350424 -0.01435817 0.00000000 -0.05024160 0.60281005 #> x7 x8 x9 x10 #> 0.53291918 0.70067227 0.69438086 0.39645481 #> #> $cure #> (Intercept) x1 x2 x3 x4 x5 #> 0.18536867 0.57611279 0.28251292 0.00000000 0.00000000 0.00000000 #> x6 x7 x8 x9 x10 #> 0.00000000 0.07310365 0.00000000 0.00000000 0.00000000 #>
coef(fit2)
#> $surv #> x1 x2 x3 x4 x5 x6 x7 #> 0.00000000 0.04206453 0.00000000 0.00000000 0.00000000 0.37514690 0.34846484 #> x8 x9 x10 #> 0.51954260 0.51481914 0.23175666 #> #> $cure #> (Intercept) x1 x2 x3 x4 x5 #> 0.16781197 0.26825151 0.08013352 0.00000000 0.00000000 0.00000000 #> x6 x7 x8 x9 x10 #> 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 #>
### regularized Cox cure model with uncertain event status =========== ## simulate a toy data set.seed(123) n_obs <- 100 p <- 10 x_mat <- matrix(rnorm(n_obs * p), nrow = n_obs, ncol = p) colnames(x_mat) <- paste0("x", seq_len(p)) surv_beta <- c(rep(0, p - 5), rep(1, 5)) cure_beta <- c(rep(1, 2), rep(0, p - 2)) dat <- simData4cure(nSubject = n_obs, lambda_censor = 0.01, max_censor = 10, survMat = x_mat, survCoef = surv_beta, cureCoef = cure_beta, b0 = 0.5, p1 = 0.95, p2 = 0.95, p3 = 0.95) ## model-fitting from given design matrices fit1 <- cox_cure_net.fit(x_mat, x_mat, dat$obs_time, dat$obs_event, surv_nlambda = 5, cure_nlambda = 5, surv_alpha = 0.8, cure_alpha = 0.8) ## model-fitting from given model formula fm <- paste(paste0("x", seq_len(p)), collapse = " + ") surv_fm <- as.formula(sprintf("~ %s", fm)) cure_fm <- surv_fm fit2 <- cox_cure_net(surv_fm, cure_fm, data = dat, time = obs_time, event = obs_event, surv_nlambda = 5, cure_nlambda = 5, surv_alpha = 0.5, cure_alpha = 0.5) ## summary of BIC's BIC(fit1)
#> df surv_df cure_df BIC #> 1 2 0 2 660.0951 #> 2 3 0 3 649.6661 #> 3 8 0 8 662.4689 #> 4 8 0 8 653.7398 #> 5 9 0 9 654.4038 #> 6 4 2 2 653.2886 #> 7 5 2 3 642.7620 #> 8 10 2 8 655.9078 #> 9 10 2 8 647.3588 #> 10 11 2 9 648.0941 #> 11 5 4 1 649.1912 #> 12 7 4 3 642.3998 #> 13 11 4 7 651.2044 #> 14 12 4 8 647.3566 #> 15 13 4 9 648.1374 #> 16 9 8 1 654.9063 #> 17 11 8 3 648.6869 #> 18 14 8 6 653.0220 #> 19 16 8 8 653.6878 #> 20 17 8 9 654.4016 #> 21 10 9 1 652.2151 #> 22 12 9 3 646.4180 #> 23 15 9 6 650.8193 #> 24 16 8 8 646.8662 #> 25 17 8 9 647.5568
BIC(fit2)
#> df surv_df cure_df BIC #> 1 2 0 2 660.1474 #> 2 3 0 3 652.6713 #> 3 7 0 7 661.2005 #> 4 8 0 8 656.0703 #> 5 8 0 8 651.2207 #> 6 4 2 2 656.1159 #> 7 5 2 3 648.6202 #> 8 10 2 8 661.5912 #> 9 10 2 8 652.5115 #> 10 11 2 9 652.3423 #> 11 4 3 1 646.9408 #> 12 6 3 3 643.6201 #> 13 11 3 8 656.7878 #> 14 11 3 8 647.7554 #> 15 12 3 9 647.6006 #> 16 9 8 1 658.2755 #> 17 11 8 3 654.7547 #> 18 15 8 7 663.5134 #> 19 16 8 8 659.1908 #> 20 16 7 9 654.4754 #> 21 10 9 1 654.1403 #> 22 12 9 3 651.1082 #> 23 16 9 7 659.9199 #> 24 17 9 8 655.5378 #> 25 18 9 9 655.3806
## list of coefficient estimates based on BIC coef(fit1)
#> $surv #> x1 x2 x3 x4 x5 x6 x7 #> 0.04679663 0.00000000 0.00000000 0.00000000 0.00000000 0.47929140 0.01572914 #> x8 x9 x10 #> 0.41052238 0.00000000 0.00000000 #> #> $cure #> (Intercept) x1 x2 x3 x4 x5 #> -0.1996682 0.3369252 0.2668648 0.0000000 0.0000000 0.0000000 #> x6 x7 x8 x9 x10 #> 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 #>
coef(fit2)
#> $surv #> x1 x2 x3 x4 x5 x6 x7 #> 0.02924518 0.00000000 0.00000000 0.00000000 0.00000000 0.40564512 0.00000000 #> x8 x9 x10 #> 0.36028994 0.00000000 0.00000000 #> #> $cure #> (Intercept) x1 x2 x3 x4 x5 #> -0.2051408 0.2676769 0.2131741 0.0000000 0.0000000 0.0000000 #> x6 x7 x8 x9 x10 #> 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 #>