-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfunction_pencal_exp.R
102 lines (90 loc) · 3.54 KB
/
function_pencal_exp.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
# ------------------------------------------------------------------------------------------------
# Contains wrapper function, utility function for running experiment using pencal library
# ------------------------------------------------------------------------------------------------
#install.packages("./package/pencal_1.2.2.tar.gz", repos = NULL, type="source")
require(tidyverse)
require(pencal)
print(paste("pencal version =", packageVersion("pencal")))
# ------------------------------------------------------------------------------------------------
# Wrapper function - Fit PRC model steps 1 to 3
# ------------------------------------------------------------------------------------------------
run_prc_steps <- function(
long.data,
surv.data,
baseline.covs,
y.names,
max.ymissing = 0.2,
n.boots = 0,
n.cores = 1,
verbose = FALSE,
penalty = "ridge",
standardize = TRUE,
pfac.base.covs = 0 # vector of binary value(s) to indicate regularization for baseline covariates
) {
# Fit PRC model steps 1 to 3
# Input:
# long.data: long data
# surv.data: survival data
# y.names: vector of candidate longitudinal variable names
# Output:
# -----------------------------------------------------------------------------------
# Step 1 of PRC-LMM: estimate the LMMs
t.start <- Sys.time()
step1 <- fit_lmms(
y.names = y.names,
fixefs = ~ age.fup, # Fixed effects
ranefs = ~ age.fup | id, # Random effects
long.data = long.data,
surv.data = surv.data,
t.from.base = Years.bl, # Name of the variable containing time from baseline in long.data
max.ymissing = max.ymissing,
n.boots = n.boots,
n.cores = n.cores)
t.end <- Sys.time()
t.step1 <- as.numeric(difftime(t.end, t.start, units = "mins"))
if(verbose){
cat("\nStep 1 completed in", t.step1, "min\n")
cat("--------------------------------------------------------------------------------\n")
}
# -----------------------------------------------------------------------------------
# Step 2 of PRC-LMM: compute the summaries
t.start <- Sys.time()
step2 <- summarize_lmms(object = step1, n.cores = n.cores)
t.end <- Sys.time()
t.step2 <- as.numeric(difftime(t.end, t.start, units = "mins"))
if(verbose){
cat("\nStep 2 completed in", t.step2, "min\n")
cat("--------------------------------------------------------------------------------\n")
}
# -----------------------------------------------------------------------------------
# Step 3 of PRC-LMM: fit the penalized Cox models
t.start <- Sys.time()
step3 <- fit_prclmm(
object = step2,
surv.data = surv.data,
baseline.covs = as.formula(paste("~", paste(baseline.covs, collapse = "+"))),
penalty = penalty,
standardize = standardize,
pfac.base.covs = pfac.base.covs, # Change it to a vector to penalize baseline covariates
n.cores = n.cores
)
t.end <- Sys.time()
t.step3 <- as.numeric(difftime(t.end, t.start, units = "mins"))
if (verbose){
cat("\nStep 3 completed in", t.step3, "min\n")
cat("--------------------------------------------------------------------------------\n")
}
# -----------------------------------------------------------------------------------
return (list(
step1 = step1,
step2 = step2,
step3 = step3,
runtimes = list(
step1 = t.step1,
step2 = t.step2,
step3 = t.step3,
total = t.step1 + t.step2 + t.step3
)
))
}
# ------------------------------------------------------------------------------------------------