R/f.sim.per.splines.R
In crtahlin/peRiodiCS: Functions for Generating Periodic Curves
#' @title Simulation function
#'
#' # TODO: change references to rcs.per to rcs_per and same for cs.per, etc
#'
#' @description Runs simulations generating many models for CS, CS.PER, RCS.PER and returns results to be able to compare them.
#'
#' @param B number of simulations
#' @param n number of points in training sample?
#' @param n.new number of points in test sample?
#' @param nk number of knots
#' @param knots vector with locations of knots
#' @param quatiles.cs quantiles at which to place the knots by default?
#' @param par1sin constant to add to "shift" the sine function to generate values in acceptable range (defaults to 1, to have function values between 0 and 2)
#' @param par2sin factor to multiply the sine function by (<1, "squashing" the function min-max difference; defaults to 0.25, to have mix&max of the output differ by 0.5 from original 2)
#' @param par3sin constant to add to additionaly "shift" the whole function to wanted level (defaults to 0.25, to raise the whole function above 0)
#' @param par4sin the number of sine cycles per one of our periods (defaults to 1 cycle)
#' @param par5sin offset of the sine curve in multiples of 2*Pi (defaults to 0); with 0.25 we would get cosine curve
#' @param tmax the end of the interval we are using, e.g. 3 for "3 years"; useful if studying trends (defaults to 1)
#' @param add_trend boolean. whether to simulate a trend in probability changeing over years. Trend also uses a sine function. defaults to FALSE.
#' @param par1trend constant to add to "shift" the trend part. (defaults to 1, to have function values between 0 and 2)
#' @param par2trend factor to multiply the trend sine function (<1, squashin it; defaults to 0.1, to have min&max differ by 0.2)
#' @param par3trend constant to add to additionaly shift the whole trend to wanted level (defaults to 0, for no additional shifting, e.g. start at 0)
#' @param par4trend the number of trend sine cycles per one of our periods (cycles) (defaults to 1/10, for repeating after 10 cycles)
#' @param par5trend offset of the trend sine curve in multiples of 2*Pi (defaults to 0, for a sine curve)
#' @param max_prob_value the maximum generated allowed probability value (defaults to 0.95)
#' @param min_prob_value the minimum allowed generated probability value (defaults to 0.05)
#' @param set_seed sets the seed - the same for each simulation "batch" (defaults to 1234)
#'
#' @export
f.sim.per.splines=function(B=100,
n=100,
n.new=100,
nk=5,
knots=NULL,
quantiles.cs = NULL, # c(0.05, 0.25, 0.50, 0.75, 0.95),
#parameters for the generation of the periodic data
par1sin=1,
par2sin=0.25,
par3sin=0.25,
par4sin=1,
par5sin=0,
# number of "cycles" to take
tmax = 1,
# parameters for trend
add_trend = FALSE,
par1trend = 1,
par2trend = 0.1,
par3trend = 0,
par4trend = 1/10,
par5trend = 0,
# min/max allowed probability
max_prob_value = 0.95,
min_prob_value = 0.05,
set_seed = 1234){
# set the seed ####
set.seed(set_seed)
############### initialization of the outputs ####
num.res = 3 + 3 + 2 + 4*2 + 2*3 + 3 + 3 # number of result columns to initialize
#4*2: AUC, 2*3: calibration, 3: LRT for training models, 3: score test for training models
my.res = matrix(NA, nrow=B, ncol=num.res)
# add names to columns
colnames(my.res) <-
c("brier.rcs.train", "brier.rcs.per.train", "brier.cs.per.train",
"brier.rcs", "brier.rcs.per", "brier.cs.per",
"AUCTrainEst.rcs", "AUCTrainEst.rcs.per", "AUCTrainEst.cs.per",
"AUCNewEst.rcs", "AUCNewEst.rcs.per", "AUCNewEst.cs.per",
"cal.rcs.1", "cal.rcs.per.1", "cal.cs.per.1",
"cal.rcs.2", "cal.rcs.per.2", "cal.cs.per.2",
"lrt.p.value.rcs.train", "lrt.p.value.rcs.per.train", "lrt.p.value.cs.per.train",
"score.p.value.rcs.train", "score.p.value.rcs.per.train", "score.p.value.cs.per.train",
"AUCTrainMax", "AUCNewMax", "prop.events.train", "prop.events.test")
# theoretical min and max of the covariate
min.th.x=0
max.th.x=1
# points where the coverage is tested ####
x.test.points = seq(min.th.x + 0.025, max.th.x - 0.025, by=.025)
num.test.points = length(x.test.points)
#matrices where there the coverage probabilities are stored
#predicted probabilities ####
prob.coverage.cs.per=prob.coverage.rcs=prob.coverage.rcs.per=matrix(NA, ncol=num.test.points*4, nrow=B)
prob.coverage.cs.per.train=prob.coverage.rcs.train=prob.coverage.rcs.per.train=matrix(NA, ncol=num.test.points*4, nrow=B)
#predicted lp ####
lp.coverage.cs.per=lp.coverage.rcs=lp.coverage.rcs.per=matrix(NA, ncol=num.test.points*4, nrow=B)
lp.coverage.cs.per.train=lp.coverage.rcs.train=lp.coverage.rcs.per.train=matrix(NA, ncol=num.test.points*4, nrow=B)
# structures to hold models
models_rcs <- list()
models_rcs.per <- list()
models_cs.per <- list()
# define probability function in one place
probability_function <-
function(x,
par1sin. = par1sin,
par2sin. = par2sin,
par3sin. = par3sin,
par4sin. = par4sin,
par5sin. = par5sin,
add_trend. = add_trend,
par1trend. = par1trend,
par2trend. = par2trend,
par3trend. = par3trend,
par4trend. = par4trend,
par5trend. = par5trend,
max_prob_value. = max_prob_value,
min_prob_value. = min_prob_value) {
sine_function(x,
par1sin=par1sin.,
par2sin=par2sin.,
par3sin=par3sin.,
par4sin=par4sin.,
par5sin=par5sin.,
add_trend = add_trend.,
par1trend = par1trend.,
par2trend = par2trend.,
par3trend = par3trend.,
par4trend = par4trend.,
par5trend = par5trend.,
max_prob_value = max_prob_value.,
min_prob_value = min_prob_value.
)
}
############### beginning of the iterations
for(b in 1:B) {
# generate x on the interval allown using uniform distribution
x = runif(n, min = 0, max = tmax) # we generate 100 x values
# simulation of the probabilities
x.transf.2 = x.transf = probability_function(x)
# but from now on, we have to transform the X to keep it inside [0,1] using modulus division, e.g. (x %% 1)
x <- x %% 1
# simulation of the binary events
y = ifelse(runif(n) < x.transf.2, 1, 0) # event happens if generated random number [0,1] smaller than simulated probability
# TODO: why not use rbinom() ?
# transformation on the linear predictor scale
x.transf.2.lp = log(x.transf.2 / (1 - x.transf.2))
# generation of the design matrices ####
# classic RCS
x.rcs = rcspline.eval(x, nk = nk, inclx = TRUE)
#saving the knots used for the splines, uses the default from rms package
my.knots = attr(x.rcs, "knots")
#periodic RCS
x.rcs.per = rcs.per(x, nk = nk, xmin = min.th.x, xmax = max.th.x)
# TODO!: check that my.knot are correct for each of the functions! and that each returns the same knot locations
# TODO!: set the default quantiles.cs to NULL? otherwise they will have knots at different locations than other variants
# DONE : set to NULL for the simulation function
# periodic cubic splines
if(is.null(quantiles.cs)){
x.cs.per = cs.per(x, nk=nk, xmin=min.th.x, xmax=max.th.x)
my.knots.cs = my.knots
} else {
my.knots.cs = as.numeric(quantile(ecdf(x), quantiles.cs))
x.cs.per = cs.per(x, knots=my.knots.cs, nk=NULL, xmin=min.th.x, xmax=max.th.x)
}
# save each of the knots
my.knots.rcs = attr(x.rcs, "knots")
my.knots.rcs.per = attr(x.rcs.per, "knots")
my.knots.cs.per = attr(x.cs.per, "knots")
# TODO: maybe later explicitly use each of these in the test set?
# TODO!: save these models to results? to plot later
# estimate models
mod.rcs.per = glm(y~x.rcs.per, family="binomial") # model binary responses with different versions of splines
mod.rcs = glm(y~x.rcs, family="binomial")
mod.cs.per = glm(y~x.cs.per, family="binomial")
# save models to results placeholder
models_rcs[[b]] <- mod.rcs
models_rcs.per[[b]] <- mod.rcs.per
models_cs.per[[b]] <- mod.cs.per
# Brier's scores
brier.rcs.per.train = mean((y - mod.rcs.per$fitted.values)^2) # Brier's score on fitted values (square of diff between real and fitted)
brier.rcs.train = mean((y-mod.rcs$fitted.values)^2)
brier.cs.per.train = mean((y-mod.cs.per$fitted.values)^2)
# calculate proportion of events happening
prop.events.train = mean(y)
################### predicted probabilities, training ####
# estimated linear predictor with se
pred.rcs.per.train = predict(mod.rcs.per, se.fit = TRUE) # on the scale of linear predictors
pred.rcs.train = predict(mod.rcs, se.fit=TRUE)
pred.cs.per.train = predict(mod.cs.per, se.fit=TRUE)
# saving the linear predictors
lp.rcs.per.train = pred.rcs.per.train$fit
lp.rcs.train = pred.rcs.train$fit
lp.cs.per.train = pred.cs.per.train$fit
# deriving the estimated probabilities
p.rcs.per.train = 1/(1+exp(-lp.rcs.per.train))
p.rcs.train = 1/(1+exp(-lp.rcs.train))
p.cs.per.train = 1/(1+exp(-lp.cs.per.train))
#estimated probabilities predictor with se
p.pred.rcs.per.train = predict(mod.rcs.per, se.fit = TRUE, type="response")
p.pred.rcs.train = predict(mod.rcs, se.fit=TRUE, type="response")
p.pred.cs.per.train = predict(mod.cs.per, se.fit=TRUE, type="response")
######## coverage #############
########### coverage predicted probabilities #######
my.index=numeric(num.test.points)
for(i in 1:num.test.points){
my.index[i]=which.min(abs(x-x.test.points[i]))}
# calculate if the real (simulated) point was covered by the confidence interval (and save some other info) for each of the variants
prob.coverage.rcs.per.train[b,]=c(p.pred.rcs.per.train$fit[my.index]-1.96*p.pred.rcs.per.train$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs.per.train$fit[my.index]+1.96*p.pred.rcs.per.train$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs.per.train$fit[my.index]-1.96*p.pred.rcs.per.train$se.fit[my.index],
p.pred.rcs.per.train$fit[my.index],
p.pred.rcs.per.train$fit[my.index]+1.96*p.pred.rcs.per.train$se.fit[my.index])
prob.coverage.rcs.train[b,]=c(p.pred.rcs.train$fit[my.index]-1.96*p.pred.rcs.train$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs.train$fit[my.index]+1.96*p.pred.rcs.train$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs.train$fit[my.index]-1.96*p.pred.rcs.train$se.fit[my.index],
p.pred.rcs.train$fit[my.index],
p.pred.rcs.train$fit[my.index]+1.96*p.pred.rcs.train$se.fit[my.index])
prob.coverage.cs.per.train[b,]=c(p.pred.cs.per.train$fit[my.index]-1.96*p.pred.cs.per.train$se.fit[my.index]<=x.transf.2[my.index] & p.pred.cs.per.train$fit[my.index]+1.96*p.pred.cs.per.train$se.fit[my.index]>=x.transf.2[my.index],
p.pred.cs.per.train$fit[my.index]-1.96*p.pred.cs.per.train$se.fit[my.index],
p.pred.cs.per.train$fit[my.index],
p.pred.cs.per.train$fit[my.index]+1.96*p.pred.cs.per.train$se.fit[my.index])
################ end predicted probabilities training
##### predicted lp coverage, training
lp.coverage.rcs.per.train[b,]=c(pred.rcs.per.train$fit[my.index]-1.96*pred.rcs.per.train$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs.per.train$fit[my.index]+1.96*pred.rcs.per.train$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs.per.train$fit[my.index]-1.96*pred.rcs.per.train$se.fit[my.index],
pred.rcs.per.train$fit[my.index],
pred.rcs.per.train$fit[my.index]+1.96*pred.rcs.per.train$se.fit[my.index])
lp.coverage.rcs.train[b,]=c(pred.rcs.train$fit[my.index]-1.96*pred.rcs.train$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs.train$fit[my.index]+1.96*pred.rcs.train$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs.train$fit[my.index]-1.96*pred.rcs.train$se.fit[my.index],
pred.rcs.train$fit[my.index],
pred.rcs.train$fit[my.index]+1.96*pred.rcs.train$se.fit[my.index])
lp.coverage.cs.per.train[b,]=c(pred.cs.per.train$fit[my.index]-1.96*pred.cs.per.train$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.cs.per.train$fit[my.index]+1.96*pred.cs.per.train$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.cs.per.train$fit[my.index]-1.96*pred.cs.per.train$se.fit[my.index],
pred.cs.per.train$fit[my.index],
pred.cs.per.train$fit[my.index]+1.96*pred.cs.per.train$se.fit[my.index])
# end predicted lp coverage, training
# AUC max on test data
AUCTrainMax = as.numeric(performance(prediction(x.transf.2, y), "auc")@y.values)
#AUC estimated on test data
AUCTrainEst.cs.per = as.numeric(performance(prediction(p.pred.cs.per.train$fit, y), "auc")@y.values)
AUCTrainEst.rcs.per=as.numeric(performance(prediction(p.pred.rcs.per.train$fit, y), "auc")@y.values)
AUCTrainEst.rcs=as.numeric(performance(prediction(p.pred.rcs.train$fit, y), "auc")@y.values)
########## lrt for variable selection
lrt.p.value.cs.per.train = 1-pchisq(-1*(mod.cs.per$deviance-mod.cs.per$null.deviance), mod.cs.per$df.null-mod.cs.per$df.residual)
lrt.p.value.rcs.per.train = 1-pchisq(-1*(mod.rcs.per$deviance-mod.rcs.per$null.deviance), mod.rcs.per$df.null-mod.rcs.per$df.residual)
lrt.p.value.rcs.train = 1-pchisq(-1*(mod.rcs$deviance-mod.rcs$null.deviance), mod.rcs$df.null-mod.rcs$df.residual)
########## score for variable selection
score.p.value.cs.per.train=anova(mod.cs.per, test="Rao")[2,6]
score.p.value.rcs.per.train=anova(mod.rcs.per, test="Rao")[2,6]
score.p.value.rcs.train=anova(mod.rcs, test="Rao")[2,6]
#4debug
x.old=x
############## new data (test set) #######################
x=runif(n.new, min = 0, max = tmax)
x.transf.2=x.transf = probability_function(x)
x <- x %% 1 # merge into one cycle
ynew=ifelse(runif(n)<x.transf, 1, 0)
x.transf.2.lp=log(x.transf.2/(1-x.transf.2))
#using the knots specified in the training data
x.rcs.per = rcs.per(x, knots=my.knots, xmin=min.th.x, xmax=max.th.x)
x.rcs = rcspline.eval(x, knots=my.knots, inclx=TRUE)
x.cs.per = cs.per(x, knots=my.knots.cs, nk=NULL, xmax=max.th.x, xmin=min.th.x)
# construct data frames to hold data to use in test models
new.data.rcs.per = data.frame(ynew=ynew, x.rcs.per)
dimnames(new.data.rcs.per)[[2]][-1]=seq(1:dim(x.rcs.per)[2])
new.data.rcs = data.frame(ynew=ynew, x.rcs)
dimnames(new.data.rcs)[[2]][-1]=seq(1:dim(x.rcs)[2])
new.data.cs.per=data.frame(ynew=ynew, x.cs.per)
dimnames(new.data.cs.per)[[2]][-1]=seq(1:dim(x.cs.per)[2])
#estimated linear predictor with se
pred.rcs.per = predict(mod.rcs.per, newdata=new.data.rcs.per, se.fit = TRUE)
pred.rcs = predict(mod.rcs, newdata=new.data.rcs, se.fit=TRUE)
pred.cs.per = predict(mod.cs.per, newdata=new.data.cs.per, se.fit=TRUE)
#saving the linear predictors
lp.rcs.per = pred.rcs.per$fit
lp.rcs = pred.rcs$fit
lp.cs.per = pred.cs.per$fit
#deriving the estimated probabilities
p.rcs.per = 1/(1+exp(-lp.rcs.per))
p.rcs = 1/(1+exp(-lp.rcs))
p.cs.per = 1/(1+exp(-lp.cs.per))
#estimated probabilities predictor with se
p.pred.rcs.per = predict(mod.rcs.per, newdata=new.data.rcs.per, se.fit = TRUE, type="response")
p.pred.rcs = predict(mod.rcs, newdata=new.data.rcs, se.fit=TRUE, type="response")
p.pred.cs.per = predict(mod.cs.per, newdata=new.data.cs.per, se.fit=TRUE, type="response")
############ brier's score on new data ###############
brier.rcs.per = mean((ynew-p.rcs.per)^2)
brier.rcs = mean((ynew-p.rcs)^2)
brier.cs.per = mean((ynew-p.cs.per)^2)
#proportion of events on new data
prop.events.test=mean(ynew)
############ Prediction CI coverage score on new data ###############
my.index=numeric(num.test.points)
for(i in 1:num.test.points){
my.index[i]=which.min(abs(x-x.test.points[i]))}
prob.coverage.rcs.per[b,]=c(p.pred.rcs.per$fit[my.index]-1.96*p.pred.rcs.per$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs.per$fit[my.index]+1.96*p.pred.rcs.per$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs.per$fit[my.index]-1.96*p.pred.rcs.per$se.fit[my.index],
p.pred.rcs.per$fit[my.index],
p.pred.rcs.per$fit[my.index]+1.96*p.pred.rcs.per$se.fit[my.index])
prob.coverage.rcs[b,]=c(p.pred.rcs$fit[my.index]-1.96*p.pred.rcs$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs$fit[my.index]+1.96*p.pred.rcs$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs$fit[my.index]-1.96*p.pred.rcs$se.fit[my.index],
p.pred.rcs$fit[my.index],
p.pred.rcs$fit[my.index]+1.96*p.pred.rcs$se.fit[my.index])
prob.coverage.cs.per[b,]=c(p.pred.cs.per$fit[my.index]-1.96*p.pred.cs.per$se.fit[my.index]<=x.transf.2[my.index] & p.pred.cs.per$fit[my.index]+1.96*p.pred.cs.per$se.fit[my.index]>=x.transf.2[my.index],
p.pred.cs.per$fit[my.index]-1.96*p.pred.cs.per$se.fit[my.index],
p.pred.cs.per$fit[my.index],
p.pred.cs.per$fit[my.index]+1.96*p.pred.cs.per$se.fit[my.index])
############# end coverage of predicted probabilities
#### coverage of linear predictor, new data
##### predicted lp coverage, training
lp.coverage.rcs.per[b,]=c(pred.rcs.per$fit[my.index]-1.96*pred.rcs.per$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs.per$fit[my.index]+1.96*pred.rcs.per$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs.per$fit[my.index]-1.96*pred.rcs.per$se.fit[my.index],
pred.rcs.per$fit[my.index],
pred.rcs.per$fit[my.index]+1.96*pred.rcs.per$se.fit[my.index])
lp.coverage.rcs[b,]=c(pred.rcs$fit[my.index]-1.96*pred.rcs$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs$fit[my.index]+1.96*pred.rcs$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs$fit[my.index]-1.96*pred.rcs$se.fit[my.index],
pred.rcs$fit[my.index],
pred.rcs$fit[my.index]+1.96*pred.rcs$se.fit[my.index])
lp.coverage.cs.per[b,]=c(pred.cs.per$fit[my.index]-1.96*pred.cs.per$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.cs.per$fit[my.index]+1.96*pred.cs.per$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.cs.per$fit[my.index]-1.96*pred.cs.per$se.fit[my.index],
pred.cs.per$fit[my.index],
pred.cs.per$fit[my.index]+1.96*pred.cs.per$se.fit[my.index])
############## end of coverage of linear predictor, new data
# TODO: are these three lines duplicated from above (not needed)?
brier.rcs.per=mean((ynew-p.rcs.per)^2)
brier.rcs=mean((ynew-p.rcs)^2)
brier.cs.per=mean((ynew-p.cs.per)^2)
#AUC max on test data
AUCNewMax = as.numeric(performance(prediction(x.transf.2, ynew), "auc")@y.values)
#AUC estimated on test data
AUCNewEst.cs.per = as.numeric(performance(prediction(p.pred.cs.per$fit, ynew), "auc")@y.values)
AUCNewEst.rcs.per = as.numeric(performance(prediction(p.pred.rcs.per$fit, ynew), "auc")@y.values)
AUCNewEst.rcs = as.numeric(performance(prediction(p.pred.rcs$fit, ynew), "auc")@y.values)
#calibration on new data
#intercept and slope, 6 new parameters
cal.cs.per=glm(ynew~pred.cs.per$fit, family="binomial")$coef
cal.rcs.per=glm(ynew~pred.rcs.per$fit, family="binomial")$coef
cal.rcs=glm(ynew~pred.rcs$fit, family="binomial")$coef
my.res[b,]=c(brier.rcs.train, brier.rcs.per.train, brier.cs.per.train,
brier.rcs, brier.rcs.per, brier.cs.per,
AUCTrainEst.rcs, AUCTrainEst.rcs.per, AUCTrainEst.cs.per,
AUCNewEst.rcs, AUCNewEst.rcs.per, AUCNewEst.cs.per,
cal.rcs[1], cal.rcs.per[1], cal.cs.per[1],
cal.rcs[2], cal.rcs.per[2], cal.cs.per[2],
lrt.p.value.rcs.train, lrt.p.value.rcs.per.train, lrt.p.value.cs.per.train,
score.p.value.rcs.train, score.p.value.rcs.per.train, score.p.value.cs.per.train,
AUCTrainMax, AUCNewMax, prop.events.train, prop.events.test)
}#end for b
# TODOs
# -pojdi skozi izvajanje, preglej rezultate
# -dodaj možnost generiranja podatkov z zamaknjeno sin funkcijo
# -dodaj možnost izvajanja z non-sin funkcijo?
# -dodaj možnost izvajanja
return(list(my.res=my.res, # results
# coverages
prob.coverage.rcs=prob.coverage.rcs,
prob.coverage.rcs.per=prob.coverage.rcs.per,
prob.coverage.cs.per=prob.coverage.cs.per,
lp.coverage.rcs=lp.coverage.rcs,
lp.coverage.rcs.per=lp.coverage.rcs.per,
lp.coverage.cs.per=lp.coverage.cs.per,
prob.coverage.rcs.train=prob.coverage.rcs.train,
prob.coverage.rcs.per.train=prob.coverage.rcs.per.train,
prob.coverage.cs.per.train=prob.coverage.cs.per.train,
lp.coverage.rcs.train=lp.coverage.rcs.train,
lp.coverage.rcs.per.train=lp.coverage.rcs.per.train,
lp.coverage.cs.per.train=lp.coverage.cs.per.train,
# points at which coverage was tested
x.test.points=x.test.points, num.test.points=num.test.points,
#saving the parameters of the simulation
B=B, n=n, n.new=n.new, nk=nk, knots=my.knots, knots.cs=my.knots.cs,
quantiles.cs=quantiles.cs,
#parameters for the generation of the periodic data
par1sin=par1sin, par2sin=par2sin, par3sin=par3sin, par4sin=par4sin, par5sin=par5sin, prop.events.train=prop.events.train, prop.events.test=prop.events.test,
x=x, x.transf.2.lp=x.transf.2.lp,
tmax = tmax,
# trend parameters
add_trend = add_trend,
par1trend = par1trend,
par2trend = par2trend,
par3trend = par3trend,
par4trend = par4trend,
par5trend = par5trend,
max_prob_value = max_prob_value,
min_prob_value = min_prob_value,
# seed
set_seed = set_seed,
# knots info explicitly
my.knots.rcs = my.knots.rcs,
my.knots.rcs.per = my.knots.rcs.per,
my.knots.cs.per = my.knots.cs.per #,
# estimated models for each simulation
# NOTE: remove temporarily, due to large size of results
# models_rcs = models_rcs,
# models_rcs.per = models_rcs.per,
# models_cs.per = models_cs.per
))
}#end function sim
crtahlin/peRiodiCS documentation built on May 15, 2019, 2:50 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' @title Simulation function
#'
#' # TODO: change references to rcs.per to rcs_per and same for cs.per, etc
#'
#' @description Runs simulations generating many models for CS, CS.PER, RCS.PER and returns results to be able to compare them.
#'
#' @param B number of simulations
#' @param n number of points in training sample?
#' @param n.new number of points in test sample?
#' @param nk number of knots
#' @param knots vector with locations of knots
#' @param quatiles.cs quantiles at which to place the knots by default?
#' @param par1sin constant to add to "shift" the sine function to generate values in acceptable range (defaults to 1, to have function values between 0 and 2)
#' @param par2sin factor to multiply the sine function by (<1, "squashing" the function min-max difference; defaults to 0.25, to have mix&max of the output differ by 0.5 from original 2)
#' @param par3sin constant to add to additionaly "shift" the whole function to wanted level (defaults to 0.25, to raise the whole function above 0)
#' @param par4sin the number of sine cycles per one of our periods (defaults to 1 cycle)
#' @param par5sin offset of the sine curve in multiples of 2*Pi (defaults to 0); with 0.25 we would get cosine curve
#' @param tmax the end of the interval we are using, e.g. 3 for "3 years"; useful if studying trends (defaults to 1)
#' @param add_trend boolean. whether to simulate a trend in probability changeing over years. Trend also uses a sine function. defaults to FALSE.
#' @param par1trend constant to add to "shift" the trend part. (defaults to 1, to have function values between 0 and 2)
#' @param par2trend factor to multiply the trend sine function (<1, squashin it; defaults to 0.1, to have min&max differ by 0.2)
#' @param par3trend constant to add to additionaly shift the whole trend to wanted level (defaults to 0, for no additional shifting, e.g. start at 0)
#' @param par4trend the number of trend sine cycles per one of our periods (cycles) (defaults to 1/10, for repeating after 10 cycles)
#' @param par5trend offset of the trend sine curve in multiples of 2*Pi (defaults to 0, for a sine curve)
#' @param max_prob_value the maximum generated allowed probability value (defaults to 0.95)
#' @param min_prob_value the minimum allowed generated probability value (defaults to 0.05)
#' @param set_seed sets the seed - the same for each simulation "batch" (defaults to 1234)
#'
#' @export
f.sim.per.splines=function(B=100,
n=100,
n.new=100,
nk=5,
knots=NULL,
quantiles.cs = NULL, # c(0.05, 0.25, 0.50, 0.75, 0.95),
#parameters for the generation of the periodic data
par1sin=1,
par2sin=0.25,
par3sin=0.25,
par4sin=1,
par5sin=0,
# number of "cycles" to take
tmax = 1,
# parameters for trend
add_trend = FALSE,
par1trend = 1,
par2trend = 0.1,
par3trend = 0,
par4trend = 1/10,
par5trend = 0,
# min/max allowed probability
max_prob_value = 0.95,
min_prob_value = 0.05,
set_seed = 1234){
# set the seed ####
set.seed(set_seed)
############### initialization of the outputs ####
num.res = 3 + 3 + 2 + 4*2 + 2*3 + 3 + 3 # number of result columns to initialize
#4*2: AUC, 2*3: calibration, 3: LRT for training models, 3: score test for training models
my.res = matrix(NA, nrow=B, ncol=num.res)
# add names to columns
colnames(my.res) <-
c("brier.rcs.train", "brier.rcs.per.train", "brier.cs.per.train",
"brier.rcs", "brier.rcs.per", "brier.cs.per",
"AUCTrainEst.rcs", "AUCTrainEst.rcs.per", "AUCTrainEst.cs.per",
"AUCNewEst.rcs", "AUCNewEst.rcs.per", "AUCNewEst.cs.per",
"cal.rcs.1", "cal.rcs.per.1", "cal.cs.per.1",
"cal.rcs.2", "cal.rcs.per.2", "cal.cs.per.2",
"lrt.p.value.rcs.train", "lrt.p.value.rcs.per.train", "lrt.p.value.cs.per.train",
"score.p.value.rcs.train", "score.p.value.rcs.per.train", "score.p.value.cs.per.train",
"AUCTrainMax", "AUCNewMax", "prop.events.train", "prop.events.test")
# theoretical min and max of the covariate
min.th.x=0
max.th.x=1
# points where the coverage is tested ####
x.test.points = seq(min.th.x + 0.025, max.th.x - 0.025, by=.025)
num.test.points = length(x.test.points)
#matrices where there the coverage probabilities are stored
#predicted probabilities ####
prob.coverage.cs.per=prob.coverage.rcs=prob.coverage.rcs.per=matrix(NA, ncol=num.test.points*4, nrow=B)
prob.coverage.cs.per.train=prob.coverage.rcs.train=prob.coverage.rcs.per.train=matrix(NA, ncol=num.test.points*4, nrow=B)
#predicted lp ####
lp.coverage.cs.per=lp.coverage.rcs=lp.coverage.rcs.per=matrix(NA, ncol=num.test.points*4, nrow=B)
lp.coverage.cs.per.train=lp.coverage.rcs.train=lp.coverage.rcs.per.train=matrix(NA, ncol=num.test.points*4, nrow=B)
# structures to hold models
models_rcs <- list()
models_rcs.per <- list()
models_cs.per <- list()
# define probability function in one place
probability_function <-
function(x,
par1sin. = par1sin,
par2sin. = par2sin,
par3sin. = par3sin,
par4sin. = par4sin,
par5sin. = par5sin,
add_trend. = add_trend,
par1trend. = par1trend,
par2trend. = par2trend,
par3trend. = par3trend,
par4trend. = par4trend,
par5trend. = par5trend,
max_prob_value. = max_prob_value,
min_prob_value. = min_prob_value) {
sine_function(x,
par1sin=par1sin.,
par2sin=par2sin.,
par3sin=par3sin.,
par4sin=par4sin.,
par5sin=par5sin.,
add_trend = add_trend.,
par1trend = par1trend.,
par2trend = par2trend.,
par3trend = par3trend.,
par4trend = par4trend.,
par5trend = par5trend.,
max_prob_value = max_prob_value.,
min_prob_value = min_prob_value.
)
}
############### beginning of the iterations
for(b in 1:B) {
# generate x on the interval allown using uniform distribution
x = runif(n, min = 0, max = tmax) # we generate 100 x values
# simulation of the probabilities
x.transf.2 = x.transf = probability_function(x)
# but from now on, we have to transform the X to keep it inside [0,1] using modulus division, e.g. (x %% 1)
x <- x %% 1
# simulation of the binary events
y = ifelse(runif(n) < x.transf.2, 1, 0) # event happens if generated random number [0,1] smaller than simulated probability
# TODO: why not use rbinom() ?
# transformation on the linear predictor scale
x.transf.2.lp = log(x.transf.2 / (1 - x.transf.2))
# generation of the design matrices ####
# classic RCS
x.rcs = rcspline.eval(x, nk = nk, inclx = TRUE)
#saving the knots used for the splines, uses the default from rms package
my.knots = attr(x.rcs, "knots")
#periodic RCS
x.rcs.per = rcs.per(x, nk = nk, xmin = min.th.x, xmax = max.th.x)
# TODO!: check that my.knot are correct for each of the functions! and that each returns the same knot locations
# TODO!: set the default quantiles.cs to NULL? otherwise they will have knots at different locations than other variants
# DONE : set to NULL for the simulation function
# periodic cubic splines
if(is.null(quantiles.cs)){
x.cs.per = cs.per(x, nk=nk, xmin=min.th.x, xmax=max.th.x)
my.knots.cs = my.knots
} else {
my.knots.cs = as.numeric(quantile(ecdf(x), quantiles.cs))
x.cs.per = cs.per(x, knots=my.knots.cs, nk=NULL, xmin=min.th.x, xmax=max.th.x)
}
# save each of the knots
my.knots.rcs = attr(x.rcs, "knots")
my.knots.rcs.per = attr(x.rcs.per, "knots")
my.knots.cs.per = attr(x.cs.per, "knots")
# TODO: maybe later explicitly use each of these in the test set?
# TODO!: save these models to results? to plot later
# estimate models
mod.rcs.per = glm(y~x.rcs.per, family="binomial") # model binary responses with different versions of splines
mod.rcs = glm(y~x.rcs, family="binomial")
mod.cs.per = glm(y~x.cs.per, family="binomial")
# save models to results placeholder
models_rcs[[b]] <- mod.rcs
models_rcs.per[[b]] <- mod.rcs.per
models_cs.per[[b]] <- mod.cs.per
# Brier's scores
brier.rcs.per.train = mean((y - mod.rcs.per$fitted.values)^2) # Brier's score on fitted values (square of diff between real and fitted)
brier.rcs.train = mean((y-mod.rcs$fitted.values)^2)
brier.cs.per.train = mean((y-mod.cs.per$fitted.values)^2)
# calculate proportion of events happening
prop.events.train = mean(y)
################### predicted probabilities, training ####
# estimated linear predictor with se
pred.rcs.per.train = predict(mod.rcs.per, se.fit = TRUE) # on the scale of linear predictors
pred.rcs.train = predict(mod.rcs, se.fit=TRUE)
pred.cs.per.train = predict(mod.cs.per, se.fit=TRUE)
# saving the linear predictors
lp.rcs.per.train = pred.rcs.per.train$fit
lp.rcs.train = pred.rcs.train$fit
lp.cs.per.train = pred.cs.per.train$fit
# deriving the estimated probabilities
p.rcs.per.train = 1/(1+exp(-lp.rcs.per.train))
p.rcs.train = 1/(1+exp(-lp.rcs.train))
p.cs.per.train = 1/(1+exp(-lp.cs.per.train))
#estimated probabilities predictor with se
p.pred.rcs.per.train = predict(mod.rcs.per, se.fit = TRUE, type="response")
p.pred.rcs.train = predict(mod.rcs, se.fit=TRUE, type="response")
p.pred.cs.per.train = predict(mod.cs.per, se.fit=TRUE, type="response")
######## coverage #############
########### coverage predicted probabilities #######
my.index=numeric(num.test.points)
for(i in 1:num.test.points){
my.index[i]=which.min(abs(x-x.test.points[i]))}
# calculate if the real (simulated) point was covered by the confidence interval (and save some other info) for each of the variants
prob.coverage.rcs.per.train[b,]=c(p.pred.rcs.per.train$fit[my.index]-1.96*p.pred.rcs.per.train$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs.per.train$fit[my.index]+1.96*p.pred.rcs.per.train$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs.per.train$fit[my.index]-1.96*p.pred.rcs.per.train$se.fit[my.index],
p.pred.rcs.per.train$fit[my.index],
p.pred.rcs.per.train$fit[my.index]+1.96*p.pred.rcs.per.train$se.fit[my.index])
prob.coverage.rcs.train[b,]=c(p.pred.rcs.train$fit[my.index]-1.96*p.pred.rcs.train$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs.train$fit[my.index]+1.96*p.pred.rcs.train$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs.train$fit[my.index]-1.96*p.pred.rcs.train$se.fit[my.index],
p.pred.rcs.train$fit[my.index],
p.pred.rcs.train$fit[my.index]+1.96*p.pred.rcs.train$se.fit[my.index])
prob.coverage.cs.per.train[b,]=c(p.pred.cs.per.train$fit[my.index]-1.96*p.pred.cs.per.train$se.fit[my.index]<=x.transf.2[my.index] & p.pred.cs.per.train$fit[my.index]+1.96*p.pred.cs.per.train$se.fit[my.index]>=x.transf.2[my.index],
p.pred.cs.per.train$fit[my.index]-1.96*p.pred.cs.per.train$se.fit[my.index],
p.pred.cs.per.train$fit[my.index],
p.pred.cs.per.train$fit[my.index]+1.96*p.pred.cs.per.train$se.fit[my.index])
################ end predicted probabilities training
##### predicted lp coverage, training
lp.coverage.rcs.per.train[b,]=c(pred.rcs.per.train$fit[my.index]-1.96*pred.rcs.per.train$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs.per.train$fit[my.index]+1.96*pred.rcs.per.train$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs.per.train$fit[my.index]-1.96*pred.rcs.per.train$se.fit[my.index],
pred.rcs.per.train$fit[my.index],
pred.rcs.per.train$fit[my.index]+1.96*pred.rcs.per.train$se.fit[my.index])
lp.coverage.rcs.train[b,]=c(pred.rcs.train$fit[my.index]-1.96*pred.rcs.train$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs.train$fit[my.index]+1.96*pred.rcs.train$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs.train$fit[my.index]-1.96*pred.rcs.train$se.fit[my.index],
pred.rcs.train$fit[my.index],
pred.rcs.train$fit[my.index]+1.96*pred.rcs.train$se.fit[my.index])
lp.coverage.cs.per.train[b,]=c(pred.cs.per.train$fit[my.index]-1.96*pred.cs.per.train$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.cs.per.train$fit[my.index]+1.96*pred.cs.per.train$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.cs.per.train$fit[my.index]-1.96*pred.cs.per.train$se.fit[my.index],
pred.cs.per.train$fit[my.index],
pred.cs.per.train$fit[my.index]+1.96*pred.cs.per.train$se.fit[my.index])
# end predicted lp coverage, training
# AUC max on test data
AUCTrainMax = as.numeric(performance(prediction(x.transf.2, y), "auc")@y.values)
#AUC estimated on test data
AUCTrainEst.cs.per = as.numeric(performance(prediction(p.pred.cs.per.train$fit, y), "auc")@y.values)
AUCTrainEst.rcs.per=as.numeric(performance(prediction(p.pred.rcs.per.train$fit, y), "auc")@y.values)
AUCTrainEst.rcs=as.numeric(performance(prediction(p.pred.rcs.train$fit, y), "auc")@y.values)
########## lrt for variable selection
lrt.p.value.cs.per.train = 1-pchisq(-1*(mod.cs.per$deviance-mod.cs.per$null.deviance), mod.cs.per$df.null-mod.cs.per$df.residual)
lrt.p.value.rcs.per.train = 1-pchisq(-1*(mod.rcs.per$deviance-mod.rcs.per$null.deviance), mod.rcs.per$df.null-mod.rcs.per$df.residual)
lrt.p.value.rcs.train = 1-pchisq(-1*(mod.rcs$deviance-mod.rcs$null.deviance), mod.rcs$df.null-mod.rcs$df.residual)
########## score for variable selection
score.p.value.cs.per.train=anova(mod.cs.per, test="Rao")[2,6]
score.p.value.rcs.per.train=anova(mod.rcs.per, test="Rao")[2,6]
score.p.value.rcs.train=anova(mod.rcs, test="Rao")[2,6]
#4debug
x.old=x
############## new data (test set) #######################
x=runif(n.new, min = 0, max = tmax)
x.transf.2=x.transf = probability_function(x)
x <- x %% 1 # merge into one cycle
ynew=ifelse(runif(n)<x.transf, 1, 0)
x.transf.2.lp=log(x.transf.2/(1-x.transf.2))
#using the knots specified in the training data
x.rcs.per = rcs.per(x, knots=my.knots, xmin=min.th.x, xmax=max.th.x)
x.rcs = rcspline.eval(x, knots=my.knots, inclx=TRUE)
x.cs.per = cs.per(x, knots=my.knots.cs, nk=NULL, xmax=max.th.x, xmin=min.th.x)
# construct data frames to hold data to use in test models
new.data.rcs.per = data.frame(ynew=ynew, x.rcs.per)
dimnames(new.data.rcs.per)[[2]][-1]=seq(1:dim(x.rcs.per)[2])
new.data.rcs = data.frame(ynew=ynew, x.rcs)
dimnames(new.data.rcs)[[2]][-1]=seq(1:dim(x.rcs)[2])
new.data.cs.per=data.frame(ynew=ynew, x.cs.per)
dimnames(new.data.cs.per)[[2]][-1]=seq(1:dim(x.cs.per)[2])
#estimated linear predictor with se
pred.rcs.per = predict(mod.rcs.per, newdata=new.data.rcs.per, se.fit = TRUE)
pred.rcs = predict(mod.rcs, newdata=new.data.rcs, se.fit=TRUE)
pred.cs.per = predict(mod.cs.per, newdata=new.data.cs.per, se.fit=TRUE)
#saving the linear predictors
lp.rcs.per = pred.rcs.per$fit
lp.rcs = pred.rcs$fit
lp.cs.per = pred.cs.per$fit
#deriving the estimated probabilities
p.rcs.per = 1/(1+exp(-lp.rcs.per))
p.rcs = 1/(1+exp(-lp.rcs))
p.cs.per = 1/(1+exp(-lp.cs.per))
#estimated probabilities predictor with se
p.pred.rcs.per = predict(mod.rcs.per, newdata=new.data.rcs.per, se.fit = TRUE, type="response")
p.pred.rcs = predict(mod.rcs, newdata=new.data.rcs, se.fit=TRUE, type="response")
p.pred.cs.per = predict(mod.cs.per, newdata=new.data.cs.per, se.fit=TRUE, type="response")
############ brier's score on new data ###############
brier.rcs.per = mean((ynew-p.rcs.per)^2)
brier.rcs = mean((ynew-p.rcs)^2)
brier.cs.per = mean((ynew-p.cs.per)^2)
#proportion of events on new data
prop.events.test=mean(ynew)
############ Prediction CI coverage score on new data ###############
my.index=numeric(num.test.points)
for(i in 1:num.test.points){
my.index[i]=which.min(abs(x-x.test.points[i]))}
prob.coverage.rcs.per[b,]=c(p.pred.rcs.per$fit[my.index]-1.96*p.pred.rcs.per$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs.per$fit[my.index]+1.96*p.pred.rcs.per$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs.per$fit[my.index]-1.96*p.pred.rcs.per$se.fit[my.index],
p.pred.rcs.per$fit[my.index],
p.pred.rcs.per$fit[my.index]+1.96*p.pred.rcs.per$se.fit[my.index])
prob.coverage.rcs[b,]=c(p.pred.rcs$fit[my.index]-1.96*p.pred.rcs$se.fit[my.index]<=x.transf.2[my.index] & p.pred.rcs$fit[my.index]+1.96*p.pred.rcs$se.fit[my.index]>=x.transf.2[my.index],
p.pred.rcs$fit[my.index]-1.96*p.pred.rcs$se.fit[my.index],
p.pred.rcs$fit[my.index],
p.pred.rcs$fit[my.index]+1.96*p.pred.rcs$se.fit[my.index])
prob.coverage.cs.per[b,]=c(p.pred.cs.per$fit[my.index]-1.96*p.pred.cs.per$se.fit[my.index]<=x.transf.2[my.index] & p.pred.cs.per$fit[my.index]+1.96*p.pred.cs.per$se.fit[my.index]>=x.transf.2[my.index],
p.pred.cs.per$fit[my.index]-1.96*p.pred.cs.per$se.fit[my.index],
p.pred.cs.per$fit[my.index],
p.pred.cs.per$fit[my.index]+1.96*p.pred.cs.per$se.fit[my.index])
############# end coverage of predicted probabilities
#### coverage of linear predictor, new data
##### predicted lp coverage, training
lp.coverage.rcs.per[b,]=c(pred.rcs.per$fit[my.index]-1.96*pred.rcs.per$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs.per$fit[my.index]+1.96*pred.rcs.per$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs.per$fit[my.index]-1.96*pred.rcs.per$se.fit[my.index],
pred.rcs.per$fit[my.index],
pred.rcs.per$fit[my.index]+1.96*pred.rcs.per$se.fit[my.index])
lp.coverage.rcs[b,]=c(pred.rcs$fit[my.index]-1.96*pred.rcs$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.rcs$fit[my.index]+1.96*pred.rcs$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.rcs$fit[my.index]-1.96*pred.rcs$se.fit[my.index],
pred.rcs$fit[my.index],
pred.rcs$fit[my.index]+1.96*pred.rcs$se.fit[my.index])
lp.coverage.cs.per[b,]=c(pred.cs.per$fit[my.index]-1.96*pred.cs.per$se.fit[my.index]<=x.transf.2.lp[my.index] & pred.cs.per$fit[my.index]+1.96*pred.cs.per$se.fit[my.index]>=x.transf.2.lp[my.index],
pred.cs.per$fit[my.index]-1.96*pred.cs.per$se.fit[my.index],
pred.cs.per$fit[my.index],
pred.cs.per$fit[my.index]+1.96*pred.cs.per$se.fit[my.index])
############## end of coverage of linear predictor, new data
# TODO: are these three lines duplicated from above (not needed)?
brier.rcs.per=mean((ynew-p.rcs.per)^2)
brier.rcs=mean((ynew-p.rcs)^2)
brier.cs.per=mean((ynew-p.cs.per)^2)
#AUC max on test data
AUCNewMax = as.numeric(performance(prediction(x.transf.2, ynew), "auc")@y.values)
#AUC estimated on test data
AUCNewEst.cs.per = as.numeric(performance(prediction(p.pred.cs.per$fit, ynew), "auc")@y.values)
AUCNewEst.rcs.per = as.numeric(performance(prediction(p.pred.rcs.per$fit, ynew), "auc")@y.values)
AUCNewEst.rcs = as.numeric(performance(prediction(p.pred.rcs$fit, ynew), "auc")@y.values)
#calibration on new data
#intercept and slope, 6 new parameters
cal.cs.per=glm(ynew~pred.cs.per$fit, family="binomial")$coef
cal.rcs.per=glm(ynew~pred.rcs.per$fit, family="binomial")$coef
cal.rcs=glm(ynew~pred.rcs$fit, family="binomial")$coef
my.res[b,]=c(brier.rcs.train, brier.rcs.per.train, brier.cs.per.train,
brier.rcs, brier.rcs.per, brier.cs.per,
AUCTrainEst.rcs, AUCTrainEst.rcs.per, AUCTrainEst.cs.per,
AUCNewEst.rcs, AUCNewEst.rcs.per, AUCNewEst.cs.per,
cal.rcs[1], cal.rcs.per[1], cal.cs.per[1],
cal.rcs[2], cal.rcs.per[2], cal.cs.per[2],
lrt.p.value.rcs.train, lrt.p.value.rcs.per.train, lrt.p.value.cs.per.train,
score.p.value.rcs.train, score.p.value.rcs.per.train, score.p.value.cs.per.train,
AUCTrainMax, AUCNewMax, prop.events.train, prop.events.test)
}#end for b
# TODOs
# -pojdi skozi izvajanje, preglej rezultate
# -dodaj možnost generiranja podatkov z zamaknjeno sin funkcijo
# -dodaj možnost izvajanja z non-sin funkcijo?
# -dodaj možnost izvajanja
return(list(my.res=my.res, # results
# coverages
prob.coverage.rcs=prob.coverage.rcs,
prob.coverage.rcs.per=prob.coverage.rcs.per,
prob.coverage.cs.per=prob.coverage.cs.per,
lp.coverage.rcs=lp.coverage.rcs,
lp.coverage.rcs.per=lp.coverage.rcs.per,
lp.coverage.cs.per=lp.coverage.cs.per,
prob.coverage.rcs.train=prob.coverage.rcs.train,
prob.coverage.rcs.per.train=prob.coverage.rcs.per.train,
prob.coverage.cs.per.train=prob.coverage.cs.per.train,
lp.coverage.rcs.train=lp.coverage.rcs.train,
lp.coverage.rcs.per.train=lp.coverage.rcs.per.train,
lp.coverage.cs.per.train=lp.coverage.cs.per.train,
# points at which coverage was tested
x.test.points=x.test.points, num.test.points=num.test.points,
#saving the parameters of the simulation
B=B, n=n, n.new=n.new, nk=nk, knots=my.knots, knots.cs=my.knots.cs,
quantiles.cs=quantiles.cs,
#parameters for the generation of the periodic data
par1sin=par1sin, par2sin=par2sin, par3sin=par3sin, par4sin=par4sin, par5sin=par5sin, prop.events.train=prop.events.train, prop.events.test=prop.events.test,
x=x, x.transf.2.lp=x.transf.2.lp,
tmax = tmax,
# trend parameters
add_trend = add_trend,
par1trend = par1trend,
par2trend = par2trend,
par3trend = par3trend,
par4trend = par4trend,
par5trend = par5trend,
max_prob_value = max_prob_value,
min_prob_value = min_prob_value,
# seed
set_seed = set_seed,
# knots info explicitly
my.knots.rcs = my.knots.rcs,
my.knots.rcs.per = my.knots.rcs.per,
my.knots.cs.per = my.knots.cs.per #,
# estimated models for each simulation
# NOTE: remove temporarily, due to large size of results
# models_rcs = models_rcs,
# models_rcs.per = models_rcs.per,
# models_cs.per = models_cs.per
))
}#end function sim
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.