testtype     = "Univ",
mtry         = mtry,
minbucket    = 0,
minsplit     = 0,
mincriterion = 0,
replace      = FALSE,
ntree        = 1000))
# specify response as nominal variable
traindata$y <- factor(traindata$y, ordered = FALSE)
# Grow RF classification (random forest consisting of classification trees)
# -------------------------------------------------------------------------
RF_classification <- cforest(y ~ .,
data     = traindata,
controls =
cforest_control(teststat     = "quad",
testtype     = "Univ",
mtry         = mtry,
minbucket    = 0,
minsplit     = 0,
mincriterion = 0,
replace      = FALSE,
ntree        = 1000))
# Compute prediction accuracy (using ranked probability score and error rate)
# ---------------------------------------------------------------------------
# obtain predicted class probabilities
pred_prob_ord <- predict(RF_ordinal, newdata = testdata, type = "prob")
pred_prob_cat <- predict(RF_classification, newdata = testdata, type = "prob")
# Ranked probability score
ncat          <- length(unique(testdata$y))
indicator_mat <- matrix(as.numeric(rep(testdata$y, each = ncat)), nrow = ncat, byrow = FALSE)
T_F           <- apply(indicator_mat, 2, function(x) as.numeric(x <= 1:ncat))
RPS_ord <- sum((sapply(pred_prob_ord, function(x) cumsum(x)) - T_F)^2)
RPS_cat <- sum((sapply(pred_prob_cat, function(x) cumsum(x)) - T_F)^2)
# obtain predicted classes
pred_ord <- as.numeric(sapply(pred_prob_ord, which.max))
pred_cat <- as.numeric(sapply(pred_prob_cat, which.max))
# Error rate
ER_ord <- mean(as.numeric(pred_ord) != as.numeric(testdata$y))
ER_cat <- mean(as.numeric(pred_cat) != as.numeric(testdata$y))
# Compute variable importance
# ---------------------------
# for RF ordinal
VIM_ER_ord  <- varimp(RF_ordinal)                      # (standard) error rate based VI
VIM_RPS_ord <- varimpRPS(RF_ordinal)                   # RPS-based VI
VIM_MSE_ord <- varimpMSE(RF_ordinal, scor = scores$y)  # MSE-based VI (scores have to be specified)
VIM_MAE_ord <- varimpMAE(RF_ordinal, scor = scores$y)  # MAE-based VI (scores have to be specified)
# for RF classification
VIM_ER_cat  <- varimp(RF_classification)                      # (standard) error rate based VI
VIM_RPS_cat <- varimpRPS(RF_classification)                   # RPS-based VI
VIM_MSE_cat <- varimpMSE(RF_classification, scor = scores$y)  # MSE-based VI (scores have to be specified)
VIM_MAE_cat <- varimpMAE(RF_classification, scor = scores$y)  # MAE-based VI (scores have to be specified)
# Return list with results on prediction accuracy and variable importance
# -----------------------------------------------------------------------
pred_imp <- list(y           = as.numeric(testdata$y),  # true response label
pred_ord    = pred_ord,                # response predicted by RF ordinal
pred_cat    = pred_cat,                # response predicted by RF classification
RPS_ord     = RPS_ord,                 # accuracy of RF ordinal in terms of RPS
RPS_cat     = RPS_cat,                 # accuracy of RF classification in terms of RPS
ER_ord      = ER_ord,                  # accuracy of RF ordinal in terms of the error rate
ER_cat      = ER_cat,                  # accuracy of RF classification in terms of the error rate
VIM_ER_ord  = VIM_ER_ord,              # error rate based VI for RF ordinal
VIM_RPS_ord = VIM_RPS_ord,             # RPS-based VI for RF ordinal
VIM_MSE_ord = VIM_MSE_ord,             # MSE-based VI for RF ordinal
VIM_MAE_ord = VIM_MAE_ord,             # MAE-based VI for RF ordinal
VIM_ER_cat  = VIM_ER_cat,              # error rate based VI for RF classification
VIM_RPS_cat = VIM_RPS_cat,             # RPS-based VI for RF classification
VIM_MSE_cat = VIM_MSE_cat,             # MSE-based VI for RF classification
VIM_MAE_cat = VIM_MAE_cat)             # MAE-based VI for RF classification
return(pred_imp)
}
perform_simulation <- function(seed,
ncat         = 6,
mixing       = 0.6,
correlations = TRUE,
nobs         = 200,
scores       = NULL,
nnoise       = 50, # must be more than 40
effect11     = 0.5,
effect12     = 0.75,
effect13     = 1,
effect2      = 1,
corstrength  = 0.8,
ncor         = 6,
mtry         = floor(sqrt(nnoise + 15))){
# use default scores if no scores are specified
if(!is.null(scores)){
if(!length(scores) == ncat) stop("Scores must have the same length as number of categories.")else scores <- list(y = scores)
}else{
scores <- list(y = 1:ncat)
}
# generate training dataset
traindata <- generate_data(seed = seed, n = nobs, ncat = ncat, mixing = mixing, correlations = correlations, nnoise = nnoise, ncor = ncor,
effect11 = effect11, effect12 = effect12, effect13 = effect13, effect2 = effect2, corstrength = corstrength)
# generate large test dataset
testdata  <- generate_data(seed = seed * 10000, n = 10000, ncat = ncat, mixing = mixing, correlations = correlations, nnoise = nnoise, ncor = ncor,
effect11 = effect11, effect12 = effect12, effect13 = effect13, effect2 = effect2, corstrength = corstrength)
# obtain and return results for generated training and test data using pre-specified scores
return(get_PA_VI(seed = seed, traindata = traindata, testdata = testdata, scores = scores, mtry = mtry))
}
setwd("Z:/tmp/Beratung/Maerte/")
data.PK<-read.csv("data_pk.csv", header=F, sep=";")
data.PK<-read.csv("data_pk_GSK.csv", header=F, sep=";")
data.PK<-data.PK[,c(1,2,3,4,5,6)]
colnames(data.PK)<-c("ID", "RTIM", "Cplasma", "RATE", "AMT", "EVID")
BREATH<-read.csv("30sLMA.csv", header=T, sep=",")
#load complete BIS data
BIS<-read.table("BISdata_GSK.csv", header=T, sep=";", dec=",")   # for this file ";" is necessary as separator
AMOUNT<-data.PK[data.PK$EVID==1,]
ID<-unique(data.PK$ID)
AMOUNT$Duration<-AMOUNT$AMT/AMOUNT$RATE
AMOUNT$Time2<-AMOUNT$RTIM+AMOUNT$Duration
# Add Time point RTIM = 0 min and RTIM 0 180 min
ZERO <- data.frame(ID = 1:20, RTIM = 0, Cplasma =0, RATE = 0, AMT = 0, EVID = 1, Duration = 0,  Time2 = 0)
T180 <- data.frame(ID = 1:20, RTIM = 180, Cplasma =0, RATE = 0, AMT = 0, EVID = 1, Duration = 0,  Time2 = 180)
AMOUNT <- rbind(AMOUNT,ZERO, T180)
AMOUNT <- AMOUNT[order(AMOUNT$ID, AMOUNT$RTIM),] #sortiert Spalten
AMOUNT$TAMT<-NA
## loop to calculate total amount of infused propofol per volunteer
for(i in ID){
subdata<-subset(AMOUNT, AMOUNT$ID == i)
for(j in 1:length(subdata$ID)){
subdata$TAMT[j]<-sum(subdata$AMT[1:j])
}
AMOUNT[AMOUNT$ID == i,]<-subdata
}
###### create subset with plasma data
PLASMA<-data.PK[data.PK$EVID==0,]
ID
data.PK$ID
head(data.PK)
setwd("Z:/cmmgrp/Silke/Bootstrap_SJ_HB_ALB/reproducible_files/R_Objects/")
load("NHANES_AIC.Rda")
names(NHANES_AIC[[1]])
varscale <- c("k = 1", "k = 4", "k = 3", "k = 4", "k = 5", "k = 1", "k = 1", "k = 1", "k = 4", "k = 1", "k = 3", "k = 4", "k = 4", "k = 4", "k = 1", "k = 1", "k = 1", "k = 1", "k = 1", "k = 1", "k = 1", "k = 1", "k = 4", "k = 1", "k = 1", "k = 1", "k = 11", "k = 1")
#graphics.off()
#pdf(file = "NHANES_AIC_ranking.pdf", height = 6, width = 12)
par(mfrow = c(1, 2), mar = c(5, 12, 0.5, 0.2))
plot(c(min(NHANES_AIC$orig_AIC), max(NHANES_AIC$orig_AIC)), c(1, length(NHANES_AIC$orig_AIC)), type="n", xlab="AIC (original sample)", yaxt="n", ylab = "")
grid(ny = c(length(NHANES_AIC$orig_AIC)+1), nx = 0)
points(sort(NHANES_AIC$orig_AIC), c(length(NHANES_AIC$orig_AIC):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(NHANES_AIC$orig_AIC)), " (", varscale[order(NHANES_AIC$orig_AIC)], ")", sep = ""), at = length(NHANES_AIC$orig_AIC):1, las = 2)
plot(c(min(rowMeans(NHANES_AIC$bootstrapped_AIC)), max(rowMeans(NHANES_AIC$bootstrapped_AIC))), c(1, length(rowMeans(NHANES_AIC$bootstrapped_AIC))),
type = "n", xlab = "Bootstrapped AIC (averaged value)", yaxt = "n", ylab = "")
grid(ny = c(length(NHANES_AIC$orig_AIC)+1), nx = 0)
points(sort(rowMeans(NHANES_AIC$bootstrapped_AIC)), c(length(rowMeans(NHANES_AIC$bootstrapped_AIC)):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(rowMeans(NHANES_AIC$bootstrapped_AIC))), " (", varscale[order(rowMeans(NHANES_AIC$bootstrapped_AIC))], ")", sep = ""), at = length(rowMeans(NHANES_AIC$bootstrapped_AIC)):1, las = 2)
#graphics.off()
#pdf(file = "NHANES_AIC_ranking.pdf", height = 6, width = 12)
par(mfrow = c(1, 2), mar = c(5, 12, 0.5, 0.2))
plot(c(min(NHANES_AIC$orig_AIC), max(NHANES_AIC$orig_AIC)), c(1, length(NHANES_AIC$orig_AIC)), type="n", xlab="AIC (original sample)", yaxt="n", ylab = "")
grid(ny = c(length(NHANES_AIC$orig_AIC)+1), nx = 0)
points(sort(NHANES_AIC$orig_AIC), c(length(NHANES_AIC$orig_AIC):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(NHANES_AIC$orig_AIC)), " (", varscale[order(NHANES_AIC$orig_AIC)], ")", sep = ""), at = length(NHANES_AIC$orig_AIC):1, las = 2)
plot(c(min(rowMeans(NHANES_AIC$bootstrapped_AIC)), max(rowMeans(NHANES_AIC$bootstrapped_AIC))), c(1, length(rowMeans(NHANES_AIC$bootstrapped_AIC))),
type = "n", xlab = "Bootstrapped AIC (averaged value)", yaxt = "n", ylab = "")
grid(ny = c(length(NHANES_AIC$orig_AIC)+1), nx = 0)
points(sort(rowMeans(NHANES_AIC$bootstrapped_AIC)), c(length(rowMeans(NHANES_AIC$bootstrapped_AIC)):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(rowMeans(NHANES_AIC$bootstrapped_AIC))), " (", varscale[order(rowMeans(NHANES_AIC$bootstrapped_AIC))], ")", sep = ""), at = length(rowMeans(NHANES_AIC$bootstrapped_AIC)):1, las = 2)
#graphics.off()
#graphics.off()
#pdf(file = "NHANES_AIC_ranking_subsample.pdf", height = 6, width = 6)
par(mfrow = c(1, 1), mar = c(5, 12, 0.5, 0.2))
plot(c(min(rowMeans(NHANES_AIC$subsampled_AIC)), max(rowMeans(NHANES_AIC$subsampled_AIC))), c(1, length(rowMeans(NHANES_AIC$subsampled_AIC))),
type = "n", xlab = "Subsampled AIC (averaged value)", yaxt = "n", ylab = "")
grid(ny = c(length(NHANES_AIC$orig_AIC)+1), nx = 0)
points(sort(rowMeans(NHANES_AIC$subsampled_AIC)), c(length(rowMeans(NHANES_AIC$subsampled_AIC)):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(rowMeans(NHANES_AIC$subsampled_AIC))), " (", varscale[order(rowMeans(NHANES_AIC$subsampled_AIC))], ")", sep = ""), at = length(rowMeans(NHANES_AIC$subsampled_AIC)):1, las = 2)
#graphics.off()
#graphics.off()
#pdf(file = "NHANES_AIC_Bias.pdf", height = 6, width = 8.5)
par(mfrow = c(1, 1), mar = c(5, 12, 0.5, 0.2))
plot(c(min(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$bootstrapped_AIC)), max(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$bootstrapped_AIC))), c(1, length(NHANES_AIC$orig_AIC)),
type = "n", xlab = "Difference in AIC (original sample) and averaged bootstrapped AIC", yaxt = "n", ylab = "")
for(i in 1:length(NHANES_AIC$orig_AIC)){
points(c(0, sort(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$bootstrapped_AIC))[i]), c(i, i), type = "l")
}
axis(2, labels = paste(names(sort(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$bootstrapped_AIC))), " (",
varscale[order(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$bootstrapped_AIC))], ")", sep = ""), at = 1:length(NHANES_AIC$orig_AIC), las = 2)
#graphics.off()
#graphics.off()
#pdf(file = "NHANES_AIC_Bias_subsample.pdf", height = 6, width = 8.5)
par(mfrow = c(1, 1), mar = c(5, 12, 0.5, 0.2))
plot(c(min(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$subsampled_AIC)), max(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$subsampled_AIC))), c(1, length(NHANES_AIC$orig_AIC)),
type = "n", xlab = "Difference in AIC (original sample) and averaged subsampled AIC", yaxt = "n", ylab = "")
for(i in 1:length(NHANES_AIC$orig_AIC)){
points(c(0, sort(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$subsampled_AIC))[i]), c(i, i), type = "l")
}
axis(2, labels = paste(names(sort(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$subsampled_AIC))), " (",
varscale[order(NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$subsampled_AIC))], ")", sep = ""), at = 1:length(NHANES_AIC$orig_AIC), las = 2)
#graphics.off()
# compute numbers in Table A3: Variable ranking for the modified NHANES data (uncomment lines to save the figure in the current working directory)
kat2 <- names(sort(NHANES_AIC$orig_AIC)[varscale[order(NHANES_AIC$orig_AIC)] == "k = 1"])
kat4 <- names(sort(NHANES_AIC$orig_AIC)[varscale[order(NHANES_AIC$orig_AIC)]  == "k = 3"])
kat5 <- names(sort(NHANES_AIC$orig_AIC)[varscale[order(NHANES_AIC$orig_AIC)]  == "k = 4"])
kat6 <- names(sort(NHANES_AIC$orig_AIC)[varscale[order(NHANES_AIC$orig_AIC)]  == "k = 5"])
kat12 <- names(sort(NHANES_AIC$orig_AIC)[varscale[order(NHANES_AIC$orig_AIC)] == "k = 11"])
# show part of table for metric and binary variables
(kat2 <- data.frame(
Original       = sapply(kat2, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)]))),
Bootstrap      = sapply(kat2, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean))))),
Bootstrap_diff = (sapply(kat2, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat2, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean)))))),
Subsample      = sapply(kat2, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean))))),
Subsample_diff = (sapply(kat2, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat2, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean)))))))
)
# show part of table for 4-category variables
(kat4 <- data.frame(
Original       = sapply(kat4, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)]))),
Bootstrap      = sapply(kat4, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean))))),
Bootstrap_diff = (sapply(kat4, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat4, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean)))))),
Subsample      = sapply(kat4, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean))))),
Subsample_diff = (sapply(kat4, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat4, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean)))))))
)
# show part of table for 4-category variables
(kat5 <- data.frame(
Original       = sapply(kat5, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)]))),
Bootstrap      = sapply(kat5, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean))))),
Bootstrap_diff = (sapply(kat5, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat5, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean)))))),
Subsample      = sapply(kat5, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean))))),
Subsample_diff = (sapply(kat5, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat5, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean)))))))
)
# show part of table for 4-category variables
(kat6 <- data.frame(
Original       = sapply(kat6, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)]))),
Bootstrap      = sapply(kat6, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean))))),
Bootstrap_diff = (sapply(kat6, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat6, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean)))))),
Subsample      = sapply(kat6, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean))))),
Subsample_diff = (sapply(kat6, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat6, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean)))))))
)
# show part of table for 4-category variables
(kat12 <- data.frame(
Original       = sapply(kat12, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)]))),
Bootstrap      = sapply(kat12, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean))))),
Bootstrap_diff = (sapply(kat12, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat12, function(z) which(z == names(sort(apply(NHANES_AIC$bootstrapped_AIC, 1, mean)))))),
Subsample      = sapply(kat12, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean))))),
Subsample_diff = (sapply(kat12, function(z) which(z == names(NHANES_AIC$orig_AIC[order(NHANES_AIC$orig_AIC)])))) - (sapply(kat12, function(z) which(z == names(sort(apply(NHANES_AIC$subsampled_AIC, 1, mean)))))))
)
par(mfrow = c(1, 2))
plot(NHANES_AIC$orig_AIC, NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$bootstrapped_AIC), xlab = "AIC (original sample)", ylab = "Difference in AIC (original - bootstrap)", main = "Bootstrap", pch = 16, cex = 0.5)
plot(NHANES_AIC$orig_AIC, NHANES_AIC$orig_AIC-rowMeans(NHANES_AIC$subsampled_AIC), xlab = "AIC (original sample)", ylab = "Difference in AIC (original - subsample)", main = "Subsample", pch = 16, cex = 0.5)
load("NHANES_original_AIC.Rda")
load("NHANES_bootstrap_AIC.Rda")
load("NHANES_subsample_AIC.Rda")
-------
# chosen number of boosting steps and variables in original sample
NHANES_original_AIC$boost.steps # = 309
length(NHANES_original_AIC$sel_coefs) - 1 # = 42 (-1 because intercept shall not be counted)
NHANES_original_AIC$boost.steps
# chosen number of boosting steps and variables in original sample
NHANES_original_AIC$boost.steps # = 309
length(NHANES_original_AIC$sel_coefs) - 1 # = 42 (-1 because intercept shall not be counted)
# mean chosen number of boosting steps in bootstrap samples
mean(sapply(NHANES_bootstrap_AIC, function(z) z$boost.steps))
# number of bootstrap samples in which a higher number of boosting steps than 309 was chosen
sum(sapply(NHANES_bootstrap_AIC, function(z) z$boost.steps > 309))
# mean number of parameters included in a model
mean(sapply(NHANES_bootstrap_AIC, function(z) length(z$sel_coefs) - 1))
# amount of bootstrap samples in which the model includes more than 42 parameters
mean(sapply(NHANES_bootstrap_AIC, function(z) (length(z$sel_coefs) - 1) > 42))
# amount of bootstrap samples in which the model includes more than 42 parameters
mean(sapply(NHANES_bootstrap_AIC, function(z) (length(z$sel_coefs) - 1) < 42))
# amount of bootstrap samples in which the model includes exactly 42 parameters
mean(sapply(NHANES_bootstrap_AIC, function(z) (length(z$sel_coefs) - 1) == 42))
# mean chosen number of boosting steps in subsamples
mean(sapply(NHANES_subsample_AIC, function(z) z$boost.steps))
# number of subsamples in which a higher number of boosting steps than 309 was chosen
sum(sapply(NHANES_subsample_AIC, function(z) z$boost.steps > 309))
# mean number of parameters included in a model
mean(sapply(NHANES_subsample_AIC, function(z) length(z$sel_coefs) - 1))
# amount of subsamples in which the model includes more than 42 parameters
mean(sapply(NHANES_subsample_AIC, function(z) (length(z$sel_coefs) - 1) > 42))
# amount of subsamples in which the model includes more than 42 parameters
mean(sapply(NHANES_subsample_AIC, function(z) (length(z$sel_coefs) - 1) < 42))
# amount of subsamples in which the model includes exactly 42 parameters
mean(sapply(NHANES_subsample_AIC, function(z) (length(z$sel_coefs) - 1) == 42))
# accuracy for models fit on bootstrap samples (in terms of MSE)
mean(sapply(NHANES_bootstrap_AIC, function(z) z$accuracy))
# accuracy for models fit on subsamples (in terms of MSE)
mean(sapply(NHANES_subsample_AIC, function(z) z$accuracy))
lablist <- c("bootstrap", "subsample")
par(mar = c(5, 4, 1, 1))
boxplot(sapply(NHANES_bootstrap_AIC, function(z) z$boost.steps), sapply(NHANES_subsample_AIC, function(z) z$boost.steps), ylab = "Optimal number of boosting steps")
abline(h = NHANES_original_AIC$boost.steps, col = "gray", lty = 2)
text(seq(0.9, 1.9, by = 1), -10, labels = lablist, srt = 45, pos = 1, xpd = TRUE, cex = 0.9)
graphics.off()
lablist <- c("bootstrap", "subsample")
par(mar = c(5, 4, 1, 1))
boxplot(sapply(NHANES_bootstrap_AIC, function(z) z$boost.steps), sapply(NHANES_subsample_AIC, function(z) z$boost.steps), ylab = "Optimal number of boosting steps")
abline(h = NHANES_original_AIC$boost.steps, col = "gray", lty = 2)
text(seq(0.9, 1.9, by = 1), -10, labels = lablist, srt = 45, pos = 1, xpd = TRUE, cex = 0.9)
#graphics.off()
#pdf(file = "NHANES_AIC_model_complexity.pdf", height = 4.5, width = 10)
par(mfrow = c(1,2))
barplot(table(sapply(NHANES_bootstrap_AIC, function(z) length(z$sel_coefs)-1))/1000, ylim = c(0, 0.12),
ylab = "Relative frequency", xlab = "Number of parameters in the boosting model", main = "Bootstrap", col = rep(c("gray90", "darkgray", "gray90"), c(9, 1, 15)))
barplot(table(sapply(NHANES_subsample_AIC, function(z) length(z$sel_coefs)-1))/1000, ylim = c(0, 0.12),
ylab = "Relative frequency", xlab = "Number of parameters in the boosting model", main = "Subsample", col = rep(c("gray90", "darkgray", "gray90"), c(19, 1, 3)))
load("NHANES_unmodified_p_values.Rda")
load("NHANES_10modified_p_values.Rda")
load("NHANES_1000modified_sig.Rda")
load("NHANES_10modified_pvalues.Rda")
# number of significant associations in the unmodified NHANES sample
sum(NHANES_unmodified_p_values$orig_pval <= 0.05) # = 17
# mean number of significant associations in unmodified NHANES bootstrap samples
mean(rowSums(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, function(z) z <= 0.05))) # = 18.4 (values always rounded to one digit)
# mean number of significant associations in unmodified NHANES subsamples
mean(rowSums(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, function(z) z <= 0.05))) # = 14.7
# number of significant associations in the modified NHANES sample
mean(sapply(1:1000, function(d) NHANES_1000modified_sig[[d]]$orig_sig)) # = 1.36
# mean number of significant associations in modified NHANES bootstrap samples
mean(sapply(1:1000, function(d)
sum(as.numeric(names(NHANES_1000modified_sig[[d]]$bootstrapped_sig)) *  as.numeric(NHANES_1000modified_sig[[d]]$bootstrapped_sig)
/sum(as.numeric(NHANES_1000modified_sig[[d]]$bootstrapped_sig))))) # = 6.12
# mean number of significant associations in modified NHANES subsamples
mean(sapply(1:1000, function(d)
sum(as.numeric(names(NHANES_1000modified_sig[[d]]$subsampled_sig)) *  as.numeric(NHANES_1000modified_sig[[d]]$subsampled_sig)
/sum(as.numeric(NHANES_1000modified_sig[[d]]$subsampled_sig))))) # = 1.40
#graphics.off()
#pdf(file = "NHANES_pvalues.pdf", height = 4, width = 7)
par(mfrow = c(1, 2))
plot(cbind(NHANES_unmodified_p_values$orig_pval, apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)), xlab = "p-value (original sample)", main = "NHANES data",
ylab = "Bootstrapped p-value (median)", ylim = c(0, 1), xlim = c(0, 1), pch = 16, cex = 0.5)
abline(c(0, 0), c(1, 1))
plot(cbind(as.numeric(sapply(1:10, function(z) NHANES_10modified_pvalues[[z]]$orig_pval)),
as.numeric(sapply(1:10, function(d) NHANES_10modified_pvalues[[d]]$bootstrapped_p_values))), xlab = "p-value (original sample)", main = "NHANES data with\n permuted response",,
ylab = "Bootstrapped p-value (median)", ylim = c(0, 1), xlim = c(0, 1), pch = 16, cex = 0.5)
abline(c(0, 0), c(1, 1))
#graphics.off()
# create Figure 3: Relative frequency of bootstrap samples with specified number of significant results when univariately testing the association between CRP level and 28 covariates
# (uncomment lines to save the figure in the current working directory)
#graphics.off()
#pdf(file = "NHANES_pvalues_no_sig.pdf", height = 4, width = 10)
par(mfrow = c(1, 2))
barplot(table(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 2, function(z)  sum(z < 0.05)))/10000,
ylab = "Relative frequency", xlab = "Number of significant associations",
ylim = c(0, 0.2), main = "NHANES data", col = rep(c("gray90", "darkgray", "gray90"), c(6, 1, 9)))
barplot(sapply(as.character(0:21), function(a)
sum(sapply(1:1000, function(z) as.numeric(NHANES_1000modified_sig[[z]]$bootstrapped_sig[a])), na.rm = TRUE))/
sum(sapply(1:1000, function(z) sum(NHANES_1000modified_sig[[z]]$bootstrapped_sig))),
ylab = "Relative frequency", xlab = "Number of significant associations",
ylim = c(0, 0.2), main = "NHANES data\n with permuted response", col = rep(c("gray90", "darkgray", "gray90"), c(1, 1, 20)))
#graphics.off()
# create Figure 4: Subsampled p-values versus original p-values (uncomment lines to save the figure in the current working directory)
#graphics.off()
#pdf(file = "NHANES_pvalues_subsample.pdf", height = 4, width = 7)
par(mfrow = c(1, 2))
plot(cbind(NHANES_unmodified_p_values$orig_pval, apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)), xlab = "p-value (original sample)", main = "NHANES data",
ylab = "Subsampled p-value (median)", ylim = c(0, 1), xlim = c(0, 1), pch = 16, cex = 0.5)
abline(c(0, 0), c(1, 1))
plot(cbind(as.numeric(sapply(1:10, function(z) NHANES_10modified_pvalues[[z]]$orig_pval)),
as.numeric(sapply(1:10, function(d) NHANES_10modified_pvalues[[d]]$subsampled_p_values))), xlab = "p-value (original sample)", main = "NHANES data with\n permuted response",,
ylab = "Subsampled p-value (median)", ylim = c(0, 1), xlim = c(0, 1), pch = 16, cex = 0.5)
abline(c(0, 0), c(1, 1))
#graphics.off()
# create Figure 5: Relative frequency of subsamples with specified number of significant results when univariately testing the association between CRP level and 28 covariates
# (uncomment lines to save the figure in the current working directory)
#graphics.off()
#pdf(file = "NHANES_pvalues_no_sig_subsample.pdf", height = 4, width = 10)
par(mfrow = c(1, 2))
barplot(table(apply(NHANES_unmodified_p_values$subsampled_p_values, 2, function(z)  sum(z < 0.05)))/10000,
ylab = "Relative frequency", xlab = "Number of significant associations",
ylim = c(0, 0.25), main = "NHANES data", col = rep(c("gray90", "darkgray", "gray90"), c(9, 1, 4)))
barplot(sapply(as.character(0:12), function(a)
sum(sapply(1:1000, function(z) as.numeric(NHANES_1000modified_sig[[z]]$subsampled_sig[a])), na.rm = TRUE))/
sum(sapply(1:1000, function(z) sum(NHANES_1000modified_sig[[z]]$subsampled_sig))),
ylim = c(0, 0.35), ylab = "Relative frequency", xlab = "Number of significant associations",
main = "NHANES data\n with permuted response", col = rep(c("gray90", "darkgray", "gray90"), c(1, 1, 11)))
#graphics.off()
# create Figure 6: Variable ranking by p-values and median bootstrapped p-value (uncomment lines to save the figure in the current working directory)
#graphics.off()
#pdf(file = "NHANES_pval_ranking.pdf", height = 6.5, width = 12)
varscale <- c("m = 2", "m = 5", "m = 4", "m = 5", "m = 6", "metric", "metric", "metric", "m = 5", "m = 2", "m = 4", "m = 5", "m = 5", "m = 5", "m = 2", "m = 2", "m = 2", "m = 2", "m = 2", "m = 2", "metric", "m = 2", "m = 5", "metric", "metric", "metric", "m = 12", "metric")
par(mfrow = c(1, 2), mar = c(5, 12, 2.5, 0.2))
plot(c(min(NHANES_unmodified_p_values$orig_pval), max(NHANES_unmodified_p_values$orig_pval)), c(1, length(NHANES_unmodified_p_values$orig_pval)), type = "n", xlab = "p-value (original sample)", yaxt = "n", ylab = "", main = "Original", xlim = c(0, 0.55))
grid(ny = c(length(NHANES_unmodified_p_values$orig_pval)+1), nx = 0)
points(sort(NHANES_unmodified_p_values$orig_pval), c(length(NHANES_unmodified_p_values$orig_pval):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(NHANES_unmodified_p_values$orig_pval)), " (",
varscale[order(NHANES_unmodified_p_values$orig_pval)], ")", sep = ""), at = length(NHANES_unmodified_p_values$orig_pval):1, las = 2)
apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)
plot(c(min(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)), max(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))), c(1, length(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))), type = "n", xlab = "Bootstrapped p-value (median)", yaxt = "n", ylab = "", main = "Bootstrap", xlim = c(0, 0.55))
grid(ny = c(length(NHANES_unmodified_p_values$orig_pval)+1), nx = 0)
points(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)), c(length(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))), " (",
varscale[order(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))], ")", sep = ""), at = length(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)):1, las = 2)
#graphics.off()
#graphics.off()
#pdf(file = "NHANES_pval_ranking_subsample.pdf", height = 6.5, width = 6)
varscale <- c("m = 2", "m = 5", "m = 4", "m = 5", "m = 6", "metric", "metric", "metric", "m = 5", "m = 2", "m = 4", "m = 5", "m = 5", "m = 5", "m = 2", "m = 2", "m = 2", "m = 2", "m = 2", "m = 2", "metric", "m = 2", "m = 5", "metric", "metric", "metric", "m = 12", "metric")
par(mfrow = c(1, 1), mar = c(5, 12, 2.5, 0.2))
apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)
plot(c(min(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)), max(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))), c(1, length(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))), type = "n", xlab = "Subsampled p-value (median)", yaxt = "n", ylab = "", main = "Subsample", xlim = c(0, 0.55))
grid(ny = c(length(NHANES_unmodified_p_values$orig_pval)+1), nx = 0)
points(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)), c(length(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)):1), cex = 0.5, col = "black", pch = 16)
axis(2, labels = paste(names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))), " (",
varscale[order(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))], ")", sep = ""), at = length(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)):1, las = 2)
#graphics.off()
# compute numbers in Table 2: Variable ranking for the unmodified NHANES data (uncomment lines to save the figure in the current working directory)
varscale <- c("m = 2", "m = 5", "m = 4", "m = 5", "m = 6", "metric", "metric", "metric", "m = 5", "m = 2", "m = 4", "m = 5", "m = 5", "m = 5", "m = 2", "m = 2", "m = 2", "m = 2", "m = 2", "m = 2", "metric", "m = 2", "m = 5", "metric", "metric", "metric", "m = 12", "metric")
kat2 <- names(sort(NHANES_unmodified_p_values$orig_pval)[varscale[order(NHANES_unmodified_p_values$orig_pval)] %in% c("m = 2", "metric")])
kat4 <- names(sort(NHANES_unmodified_p_values$orig_pval)[varscale[order(NHANES_unmodified_p_values$orig_pval)]  == "m = 4"])
kat5 <- names(sort(NHANES_unmodified_p_values$orig_pval)[varscale[order(NHANES_unmodified_p_values$orig_pval)]  == "m = 5"])
kat6 <- names(sort(NHANES_unmodified_p_values$orig_pval)[varscale[order(NHANES_unmodified_p_values$orig_pval)]  == "m = 6"])
kat12 <- names(sort(NHANES_unmodified_p_values$orig_pval)[varscale[order(NHANES_unmodified_p_values$orig_pval)] == "m = 12"])
# show part of table for metric and binary variables
(kat2 <- data.frame(
Original       = sapply(kat2, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)]))),
Bootstrap      = sapply(kat2, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))))),
Bootstrap_diff = (sapply(kat2, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat2, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)))))),
Subsample      = sapply(kat2, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))))),
Subsample_diff = (sapply(kat2, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat2, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)))))))
)
# show part of table for 4-category variables
(kat4 <- data.frame(
Original       = sapply(kat4, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)]))),
Bootstrap      = sapply(kat4, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))))),
Bootstrap_diff = (sapply(kat4, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat4, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)))))),
Subsample      = sapply(kat4, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))))),
Subsample_diff = (sapply(kat4, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat4, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)))))))
)
# show part of table for 5-category variables
(kat5 <- data.frame(
Original       = sapply(kat5, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)]))),
Bootstrap      = sapply(kat5, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))))),
Bootstrap_diff = (sapply(kat5, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat5, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)))))),
Subsample      = sapply(kat5, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))))),
Subsample_diff = (sapply(kat5, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat5, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)))))))
)
# show part of table for 6-category variables
(kat6 <- data.frame(
Original       = sapply(kat6, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)]))),
Bootstrap      = sapply(kat6, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))))),
Bootstrap_diff = (sapply(kat6, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat6, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)))))),
Subsample      = sapply(kat6, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))))),
Subsample_diff = (sapply(kat6, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat6, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)))))))
)
# show part of table for 12-category variables
(kat12 <- data.frame(
Original       = sapply(kat12, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)]))),
Bootstrap      = sapply(kat12, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median))))),
Bootstrap_diff = (sapply(kat12, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat12, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$bootstrapped_p_values, 1, median)))))),
Subsample      = sapply(kat12, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median))))),
Subsample_diff = (sapply(kat12, function(z) which(z == names(NHANES_unmodified_p_values$orig_pval[order(NHANES_unmodified_p_values$orig_pval)])))) - (sapply(kat12, function(z) which(z == names(sort(apply(NHANES_unmodified_p_values$subsampled_p_values, 1, median)))))))
)
#library(xtable)
#xtable(rbind(kat2, kat4, kat5, kat6, kat12))
# compute numbers in Table A2: Variable ranking for the first modified NHANES data (uncomment lines to save the figure in the current working directory)
i <- 1 # set number to 1, 2, ..., or 10 to report results for the 10th permuted datasets (in the paper we show results for the first one as an example)
kat2 <- rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)][varscale[order(NHANES_10modified_pvalues[[i]]$orig_pval)] %in% c("m = 2", "metric")]
kat4 <- rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)][varscale[order(NHANES_10modified_pvalues[[i]]$orig_pval)] == "m = 4"]
kat5 <- rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)][varscale[order(NHANES_10modified_pvalues[[i]]$orig_pval)] == "m = 5"]
kat6 <- rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)][varscale[order(NHANES_10modified_pvalues[[i]]$orig_pval)] == "m = 6"]
kat12 <- rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)][varscale[order(NHANES_10modified_pvalues[[i]]$orig_pval)] == "m = 12"]
# show part of table for metric and binary variables
(kat2 <- data.frame(
Original       = sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)])),
Bootstrap      = sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)])),
Bootstrap_diff = (sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)]))),
Subsample      = sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)])),
Subsample_diff = (sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat2, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)]))))
)
# show part of table for 4-category variables
(kat4 <- data.frame(
Original       = sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)])),
Bootstrap      = sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)])),
Bootstrap_diff = (sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)]))),
Subsample      = sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)])),
Subsample_diff = (sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat4, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)]))))
)
# show part of table for 5-category variables
(kat5 <- data.frame(
Original       = sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)])),
Bootstrap      = sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)])),
Bootstrap_diff = (sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)]))),
Subsample      = sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)])),
Subsample_diff = (sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat5, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)]))))
)
# show part of table for 6-category variables
(kat6 <- data.frame(
Original       = sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)])),
Bootstrap      = sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)])),
Bootstrap_diff = (sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)]))),
Subsample      = sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)])),
Subsample_diff = (sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat6, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)]))))
)
# show part of table for 12-category variables
(kat12 <- data.frame(
Original       = sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)])),
Bootstrap      = sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)])),
Bootstrap_diff = (sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$bootstrapped_p_values)]))),
Subsample      = sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)])),
Subsample_diff = (sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$orig_pval)]))) - (sapply(kat12, function(z) which(z == rownames(NHANES_10modified_pvalues[[i]])[order(NHANES_10modified_pvalues[[i]]$subsampled_p_values)]))))
)
library(party)
?cforest
library(TH.data)
data("mammoexp", package = "TH.data")
\end{document}
install.packages(knitr)
install.packages("knitr")
library(knitr)
library(TH.data)
data("mammoexp", package = "TH.data")
setwd("Z:/cmmgrp/Silke/RFordinal_SJ_GT_ALB/reproducible_files/")
source('novel_VIMs.R')
#source('novel_VIMs.R')
method(party)
methods(party)
?mammoexp
levels(mammoexp$ME)
summary(mammoexp$ME)
is.ordered(mammoexp$ME)
?is.ordered
?factor
