andreas/PAaSO
Watch
1
0
Fork
You've already forked PAaSO
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

transfer from old repo

Browse source
This commit is contained in:
Andreas Gammelgaard Damsbo 2026-08-19 09:27:27 +02:00
commit 277e2b8cf3
No known key found for this signature in database
1111 changed files with 83736 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

251
1 PA Decline/Til DDV/00master.R Normal file
View file

@ -0,0 +1,251 @@
##
## Master script
##
## Based on the assignment work from the ISL-course
##
## Generation 2 - 02.december.2022
## Code preparation for analysis on Denmarks Statistics server with enriched data set.
##
## Analysis plan:
## Table 1
## Figure 1: Sankey plot (drop & hop colored)
## Table 2: Linear regression model of pase_6~.
## Table 3: Elastic net prediction models of drop and hop. Performance measures referenced in text.
##
## A Rmarkdown file could be created to write the initial report with main results.
## This code is a bit of a mess, as it is the result of several iterations. It works however.
##
## ====================================================================
# Step 0: Primary outcome
## ====================================================================
# Script to run as hop and drop
pout <- "drop" # Drop to first quartile
# decl_rel
# decl_abs
# drop
# hop
## ====================================================================
## Data
## ====================================================================
setwd("/Users/au301842/PhysicalActivityandStrokeOutcome/1 PA Decline/")
source("data_set.R")
source("data_format.R")
## ====================================================================
##
## Baseline - by PASE group
##
## ====================================================================
ts_q <- X_tbl |>
select(vars) |>
mutate(pase_0_cut = factor(quantile_cut(pase_0, groups = 4)[[1]],ordered = TRUE)) |>
select(-pase_6,-pase_0) |>
tbl_summary(missing = "no",
by="pase_0_cut",
value = list(where(is.factor) ~ "2"),
type = list(mrs_0 ~ "categorical",
all_continuous() ~ "continuous2"),
statistic = list(all_continuous() ~ c("{N_nonmiss}",
"{median} ({p25}, {p75})",
"{min}, {max}",
"{mean} ({sd})"))
) |>
add_overall() |>
add_n ()
ts_q
tbl_one_rtf <- file("table1.RTF", "w")
writeLines(ts_q%>%as_gt()%>%as_rtf(), tbl_one_rtf)
close(tbl_one_rtf)
## ====================================================================
# Drops and hops
## ====================================================================
# TRUEs are patients dropping
table(X_tbl$pase_0_cut!="1"&X_tbl$pase_6_cut=="1")/nrow(X_tbl[X_tbl$pase_0_cut!="1",])
# TRUEs are percentage of patients inactive before stroke being more active after
table(X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut!="1")/nrow(X_tbl[X_tbl$pase_0_cut=="1",])
# TRUEs are percentage of patients being more active after that were inactive before stroke
table(X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut!="1")/nrow(X_tbl[X_tbl$pase_6_cut!="1",])
# Difference between hop/no-hop
t.test(X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut!="1","pase_0"],X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_0"])
summary(X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut!="1","pase_0"])
summary(X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_0"])
boxplot(X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut!="1","pase_0"],X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_0"])
# Stationary low
t.test(X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_0"],X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_6"])
boxplot(X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_0"],X_tbl[X_tbl$pase_0_cut=="1"&X_tbl$pase_6_cut=="1","pase_6"])
## ====================================================================
# Sankey plot
## ====================================================================
# source("sankey.R")
# p_delta
## ====================================================================
# Six months PASE: Bivariate and multivariate analyses
## ====================================================================
dta_lmreg <- X_tbl |>
select(vars) |>
mutate(mrs_0=factor(ifelse(mrs_0==1,1,2)))
Hmisc::label(dta_lmreg$mrs_0) <- "Pre-stroke mRS >0"
uv_reg <- tbl_uvregression(data=dta_lmreg,
method=lm,
y="pase_6",
show_single_row = where(is.factor),
estimate_fun = ~style_sigfig(.x,digits = 3),
pvalue_fun = ~style_pvalue(.x, digits = 3)
)
mu_reg <- dta_lmreg |>
lm(formula=pase_6~.,data=_) |>
tbl_regression(show_single_row = where(is.factor),
estimate_fun = ~style_sigfig(.x,digits = 3),
pvalue_fun = ~style_pvalue(.x, digits = 3)
)|>
add_n()
tbl_merge(list(uv_reg,mu_reg))
## ====================================================================
##
## Data variance
##
## Illustrating principal components.
##
## ====================================================================
# source("PCA.R")
#
#
# pca22
# ggsave("pc_plot.png",width = 18, height = 12, dpi = 300, limitsize = TRUE, units = "cm")
## ====================================================================
##
## Models
##
## ====================================================================
# source("assign_full.R")
ls <- list()
for (i in c("drop","hop")){
pout <- i
source("data_format.R")
source("regularisation_steps.R")
}
# Loop to run regularised model on both drop and hop.
# Saved in list for printing and exporting the plot.
## ====================================================================
# Step 1: data merge
## ====================================================================
tbl<-merge(ls$drop$RegularisedCoefs$'_data',ls$hop$RegularisedCoefs$'_data',by="name",all.x=T, sort=F)
## ====================================================================
# Step 2: table
## ====================================================================
com_coef_tbl<-tbl%>%
gt()%>%
fmt_number(
columns=colnames(tbl)[sapply(tbl,is.numeric)], ## Selecting all numeric
rows = everything(),
decimals = 3)%>%
tab_spanner(
label = "DROP",
columns = 2:5
)%>%
tab_spanner(
label = "HOP",
columns = 6:9
)%>%
tab_header(
title = "Model coefficients",
subtitle = "Combined table of both full and regularised model coefficients"
)
# paste0("Regularised model, (a=",
# best_alph,
# ", l=",
# round(best_lamb,3),
# ")")
com_coef_tbl
## ====================================================================
# Step 3: export
## ====================================================================
com_coef_rtf <- file("table2.RTF", "w")
writeLines(com_coef_tbl%>%as_rtf(), com_coef_rtf)
close(com_coef_rtf)
## ====================================================================
##
## Model performance
##
## Table with performance meassures for the two different models.
##
## ====================================================================
## ====================================================================
# Step 1: data set
## ====================================================================
tbl<-data.frame(Meassure=c(names(ls$drop$ConfusionMatrx$byClass),"Mean AUC"),
"Drop"=round(c(ls$drop$ConfusionMatrx$byClass,ls$drop$AUROC["Mean"]),3),
"Hop"=round(c(ls$hop$ConfusionMatrx$byClass,ls$hop$AUROC["Mean"]),3))
## ====================================================================
# Step 2: table
## ====================================================================
tbl_perf<-tbl%>%
gt()%>%
tab_header(
title = "Performance meassures",
subtitle = "Combined table of both drop and hop"
)
tbl_perf
## ====================================================================
# Step 3: export
## ====================================================================
tbl_perf_rtf <- file("table3.RTF", "w")
writeLines(tbl_perf%>%as_rtf(), tbl_perf_rtf)
close(tbl_perf_rtf)

50
1 PA Decline/Til DDV/data_format.R Normal file
View file

@ -0,0 +1,50 @@
## Article 1 outcome group definition script
## To be enriched from Statistics Denmark
##
## Based on the ItMLiHSmar2022 course
library(Hmisc)
library(dplyr)
library(daDoctoR)
library(tidyselect)
# Setting final primary output from "pout"
if (pout=="drop"){
X_tbl <- X_tbl|>
mutate(group=pase_drop_fac)
# print(quantile(as.numeric(X_tbl$pase_0)))
# print(quantile(as.numeric(X_tbl$pase_6)))
# print(summary(X_tbl$pase_0_cut))
X_tbl_f <- X_tbl|>
filter(pase_0_cut!=1)|>
select(-starts_with("pase_"))
}
if (pout=="hop"){
X_tbl <- X_tbl|>
mutate(group=pase_hop_fac)
# print(quantile(as.numeric(X_tbl$pase_0)))
# print(quantile(as.numeric(X_tbl$pase_6)))
# print(summary(X_tbl$pase_0_cut))
X_tbl_f <- X_tbl|>
filter(pase_6_cut!=1)|>
select(-starts_with("pase_"))
}
# Dropping non-complete for analysis
Xy <- X_tbl_f|>
na.omit()|> # Keeping only complete observations
select(-c(tci) # Left out of model as no present in drop-group
)|>
mutate(mrs_0=factor(ifelse(mrs_0==1,1,2))) # Sets binary mRS 0 to include in glmnet, 0 or above
label(Xy) = as.list(var.labels[match(names(Xy), names(var.labels))])
X<-dplyr::select(Xy,-c(group, -starts_with("pase_")) # Exclude primary outcome
)
y<-Xy$group

193
1 PA Decline/Til DDV/data_set.R Normal file
View file

@ -0,0 +1,193 @@
## Article 1 data set definition
## To be enriched from Statistics Denmark
##
## Based on the ItMLiHSmar2022 course
library(Hmisc)
library(dplyr)
library(daDoctoR)
library(tidyverse)
library(patchwork)
library(caret)
library(glmnet)
library(leaps)
library(pROC)
library(gt)
library(gtsummary)
library(glue)
# library(ggdendro)
library(corrplot)
## ====================================================================
# Step 2: Selection
## ====================================================================
export<-export[,c("pase_0",
"age",
"sex",
"civil",
"smoke_ever",
"smoker",
"rtreat",
"alc",
"afli",
"hypertension",
"diabetes",
"mrs_0",
"nihss_c",
"thrombolysis",
"pad",
"thrombechtomy",
"ami",
"tci",
"pase_6")]
## ====================================================================
# Step 3: Formatting variables
## ====================================================================
export$diabetes[is.na(export$diabetes)]<-"no"
export$diabetes[is.na(export$hypertension)]<-"no"
export$thrombolysis[is.na(export$thrombolysis)]<-"no"
export$thrombechtomy[is.na(export$thrombechtomy)]<-"no"
export$pad[is.na(export$pad)]<-"no"
export$ami[is.na(export$ami)]<-"no"
# export$smoker_prev <- ifelse(export$smoker=="3","yes","no")
export$smoker <- ifelse(export$smoker=="1","yes","no")
export$smoker[is.na(export$smoker)] <- "no"
# export$mrs_0[export$mrs_0==3]<-NA
dta <- export %>%
# as_tibble()%>%
mutate(any_rep=factor(ifelse(thrombolysis=="yes"|thrombechtomy=="yes","yes","no")), # If not noted, no therapy was received
male_sex= factor(ifelse(sex=="female","no","yes")),
# smoke_ever=factor(ifelse(smoke_ever=="never","no","yes")),
civil=factor(ifelse(civil=="partner","no","yes")), # Sets "yes" for not-cohabiting
rtreat=factor(ifelse(rtreat=="Placebo","no","yes")), # "Yes" receives active treatment
alc=factor(ifelse(alc=="more","yes","no")), # Yes for more than guideline
pase_0=as.numeric(pase_0),
pase_6=as.numeric(pase_6),
across(c("diabetes",
"hypertension",
"smoker",
"afli",
"pad",
"ami",
"tci",
"mrs_0"),as.factor),
across(c("nihss_c",
"age"),as.numeric )
)%>%
select(-c(sex))
## ====================================================================
# Step 4: Defining outcome
## ====================================================================
## Changed to step 7
## This is to perform proper quantile split based on actually included.
## ====================================================================
# Step 5: Ordering variables
## ====================================================================
vars <- c("age",
"male_sex",
"civil",
"pase_0",
"smoker",
"alc",
"afli",
"hypertension",
"diabetes",
"pad",
"ami",
"tci",
"mrs_0",
"nihss_c",
"any_rep",
"rtreat",
"pase_6")
dta<-dta[vars]
## ====================================================================
# Step 6: Labeling
## ====================================================================
var.labels = c(age="Age",
male_sex="Male",
civil="Living alone",
pase_0="Pre-stroke PASE score",
pase_6="Six month PASE score",
smoker="Daily or occasinally smoking",
alc="More alcohol than recommendation",
afli="AFIB",
hypertension="Hypertension",
diabetes="Diabetes",
pad="PAD",
ami="Previous MI",
tci="Previous TIA",
mrs_0="Pre-stroke mRS [-1]",
nihss_c="Acute NIHSS score",
thrombolysis="Acute thrombolysis",
thrombechtomy="Acute thrombechtomy",
any_rep="Any reperfusion therapy",
rtreat="Active trial treatment",
pase_drop_fac="PASE first quartile drop F",
pase_hop_fac="PASE first quartile hop F",
pase_0_cut="PASE 0 quartiles",
pase_6_cut="PASE 6 quartiles")
## ====================================================================
# Step 7: final data export
## ====================================================================
data_summary<-summary(dta)
# Saving "old" factorised variables
sel<-sapply(dta,is.factor)
# Reformatting factors as 1/2 for analysis
dta<-dta |>
mutate(across(where(is.factor), as.numeric))|> # Turning factors into 1(no) or 2(yes) for model. Numbered alphabetically.
mutate(across(matches(colnames(dta)[sel]), as.factor),
across(starts_with("pase_"), as.numeric))
# Filtering out non-PASE
X_tbl<-dta |>
filter(!is.na(pase_0),!is.na(pase_6))
nrow(X_tbl)
# Defining possible outcome meassures. Keeping in df for characterisation
X_tbl <- X_tbl|>
mutate(## Relative decline
pase_diff=(pase_0-pase_6),
pase_decl_rel = pase_diff/pase_0*100,
# pase_decl_rel_fac=factor(ifelse(pase_decl_rel>=rel_dif,"yes","no")),
## Absolute decline
# pase_decl_abs_fac=factor(ifelse(pase_diff>=abs_dif,"yes","no")),
## Drop
pase_0_cut=quantile_cut(as.numeric(pase_0),
groups=4,
group.names = c(as.character(1:4)),
y=as.numeric(pase_0),
ordered.f = TRUE,
inc.outs = TRUE,
detail.lst=FALSE),
pase_6_cut=quantile_cut(as.numeric(pase_6),
groups=4,
group.names = c(as.character(1:4)),
y=as.numeric(pase_0),
ordered.f = TRUE,
inc.outs = TRUE,
detail.lst=FALSE),
pase_drop_fac=factor(ifelse(pase_6_cut==1&pase_0_cut!=1,"yes","no")),
pase_hop_fac=factor(ifelse(pase_6_cut!=1&pase_0_cut==1,"yes","no")))
Hmisc::label(X_tbl) = as.list(var.labels[match(names(X_tbl), names(var.labels))])

117
1 PA Decline/Til DDV/regular_fun.R Normal file
View file

@ -0,0 +1,117 @@
## ItMLiHSmar2022
## regular_fun.R, child script
## Regularisation model building function
## Andreas Gammelgaard Damsbo, agdamsbo@clin.au.dk
##
## Now modified to use in publication
##
regular_fun<-function(X,y,K,lambdas,alpha){
n<-nrow(X)
set.seed(321)
# Using caret function to ensure both levels represented in all folds
c<-createFolds(y=y, k = K, list = FALSE, returnTrain = TRUE)
B<-yhatTestProbKeep<-list()
accTrain<-accTest<-err_train<-err_test<-auc_train<-auc_test<-matrix(nrow = K,ncol = length(lambdas))
catinfo<-levels(y)
cMatTrain<-cMatTest<-table(true=factor(c(0,0),levels=catinfo),pred=factor(c(0,0),levels=catinfo))
## Iterate over partitions
for (idx1 in 1:K){
# Status
cat('Processing fold', idx1, 'of', K,'\n')
# idx1=1
# Get training- and test sets
I_train = c!=idx1 ## Creating selection vector of TRUE/FALSE
I_test = !I_train
Xtrain = X[I_train,]
ytrain = y[I_train]
Xtest = X[I_test,]
ytest = y[I_test]
## Model matrices for glmnet
## Using the complicated approach not to include first level.
# Xmat.train<-model.matrix(~ .-1, data=Xtrain,
# contrasts.arg = lapply(Xtrain[,sapply(Xtrain, is.factor)],
# contrasts, contrasts=T))
# Xmat.test<-model.matrix(~ .-1, data=Xtest,
# contrasts.arg = lapply(Xtest[,sapply(Xtest, is.factor)],
# contrasts, contrasts=T))
# Xmat.train<-model.matrix(~.-1,Xtrain)
# Xmat.test<-model.matrix(~.-1,Xtest)
# Weights
ytrain_weight<-as.vector(1 - (table(ytrain)[ytrain] / length(ytrain)))
# ytest_weight<-as.vector(1 / (table(ytest)[ytest] / length(ytest)))
# Fit regularized linear regression model
mod<-glmnet(Xtrain, ytrain,
alpha = alpha, ## Alpha = 1 for lasso
lambda = lambdas, ## Setting lambdas
standardize = TRUE, ## Scales and centers
weights = ytrain_weight,
family = "binomial"
)
# Keep coefficients for plot
B[[idx1]] <- as.matrix(coef(mod))
# Iterate over regularization strengths to compute training- and test
# errors for individual regularization strengths.
for (idx2 in 1:length(lambdas)){
# idx2=1
# Predict
yhatTrainProb<-predict(mod,
s = lambdas[idx2],
newx = data.matrix(Xtrain),
type = "response"
)
yhatTestProb<-predict(mod,
s = lambdas[idx2],
newx = data.matrix(Xtest),
type = "response"
)
# Compute training and test error
yhatTrain = round(yhatTrainProb)
yhatTest = round(yhatTestProb)
# Make predictions categorical again (instead of 0/1 coding)
yhatTrainCat = factor(round(yhatTrainProb),levels=c("0","1"),labels=catinfo,ordered = TRUE)
yhatTestCat = factor(round(yhatTestProb),levels=c("0","1"),labels=catinfo,ordered = TRUE)
# Evaluate classifier performance
# Accuracy
# accTrain[idx1,idx2] <- sum(yhatTrainCat==ytrain)/length(ytrain)
# accTest [idx1,idx2] <- sum(yhatTestCat==ytest)/length(ytest)
# #
# # Error rate
# err_train[idx1,idx2] = 1 - accTrain[idx1,idx2]
# err_test [idx1,idx2] = 1 - accTest[idx1,idx2]
# AUROC
suppressMessages(
auc_train[idx1,idx2]<-auc(ytrain, yhatTrainCat))
suppressMessages(
auc_test [idx1,idx2]<-auc(ytest, yhatTestCat))
# Compute confusion matrices
cMatTrain = cMatTrain + table(true=ytrain,pred=yhatTrainCat)
cMatTest = cMatTest + table(true=ytest,pred=yhatTestCat)
}
}
ls<-list(mod=mod,B=B,auc_train=auc_train,auc_test=auc_test,cMatTrain=cMatTrain,cMatTest=cMatTest)
return(ls)
}

150
1 PA Decline/Til DDV/regularisation_steps.R Normal file
View file

@ -0,0 +1,150 @@
## ItMLiHSmar2022
## regularisation_steps.R, child script
## Regularised model building and analysation for assignment
## Andreas Gammelgaard Damsbo, agdamsbo@clin.au.dk
##
## Now modified to use in publication
##
## ====================================================================
## Step 0: data import and wrangling
## ====================================================================
setwd("/Users/au301842/PhysicalActivityandStrokeOutcome/1 PA Decline/")
# source("data_format.R")
y1<-factor(as.integer(y)-1) ## Outcome is required to be factor of 0 or 1.
## ====================================================================
## Step 1: settings
## ====================================================================
## Folds
K=10
set.seed(3)
c<-caret::createFolds(y=y,
k = K,
list = FALSE,
returnTrain = TRUE) # Foldids for alpha tuning
## Defining tuning parameters
lambdas=2^seq(-10, 5, 1)
alphas<-seq(0,1,.1)
## Weights for models
weighted=TRUE
if (weighted == TRUE) {
wght<-as.vector(1 - (table(y)[y] / length(y)))
} else {
wght <- rep(1, nrow(y))
}
## Standardise numeric
## Centered and
## ====================================================================
## Step 2: all cross validations for each alpha
## ====================================================================
library(furrr)
library(purrr)
library(doMC)
registerDoMC(cores=6)
# Nested CVs with analysis for all lambdas for each alpha
#
set.seed(3)
cvs <- future_map(alphas, function(a){
cv.glmnet(model.matrix(~.-1,X),
y1,
weights = wght,
lambda=lambdas,
type.measure = "deviance", # This is standard measure and recommended for tuning
foldid = c, # Per recommendation the folds are kept for alpha optimisation
alpha=a,
standardize=TRUE,
family=quasibinomial,
keep=TRUE) # Same as binomial, but not as picky
})
## ====================================================================
# Step 3: optimum lambda for each alpha
## ====================================================================
# For each alpha, lambda is chosen for the lowest meassure (deviance)
each_alpha <- sapply(seq_along(alphas), function(id) {
each_cv <- cvs[[id]]
alpha_val <- alphas[id]
index_lmin <- match(each_cv$lambda.min,
each_cv$lambda)
c(lamb = each_cv$lambda.min,
alph = alpha_val,
cvm = each_cv$cvm[index_lmin])
})
# Best lambda
best_lamb <- min(each_alpha["lamb", ])
# Alpha is chosen for best lambda with lowest model deviance, each_alpha["cvm",]
best_alph <- each_alpha["alph",][each_alpha["cvm",]==min(each_alpha["cvm",]
[each_alpha["lamb",] %in% best_lamb])]
## https://stackoverflow.com/questions/42007313/plot-an-roc-curve-in-r-with-ggplot2
p_roc<-roc.glmnet(cvs[[1]]$fit.preval, newy = y)[[match(best_alph,alphas)]]|> # Plots performance from model with best alpha
ggplot(aes(FPR,TPR)) +
geom_step() +
coord_cartesian(xlim=c(0,1), ylim=c(0,1)) +
geom_abline()+
theme_bw()
## ====================================================================
# Step 4: Creating the final model
## ====================================================================
source("regular_fun.R") # Custom function
optimised_model<-regular_fun(X,y1,K,lambdas=best_lamb,alpha=best_alph)
# With lambda and alpha specified, the function is just a k-fold cross-validation wrapper,
# but keeps model performance figures from each fold.
list2env(optimised_model,.GlobalEnv)
# Function outputs a list, which is unwrapped to Env.
# See source script for reference.
## ====================================================================
# Step 5: creating table of coefficients for inference
## ====================================================================
Bmatrix<-matrix(unlist(B),ncol=10)
Bmedian<-apply(Bmatrix,1,median)
Bmean<-apply(Bmatrix,1,mean)
reg_coef_tbl<-tibble(
name = c("Intercept",Hmisc::label(X)),
medianX = round(Bmedian,5),
ORmed = round(exp(Bmedian),5),
meanX = round(Bmean,5),
ORmea = round(exp(Bmean),5))%>%
# arrange(desc(abs(medianX)))%>%
gt()
## ====================================================================
# Step 6: plotting predictive performance
## ====================================================================
reg_cfm<-confusionMatrix(cMatTest)
reg_auc_sum<-summary(auc_test[,1])
## ====================================================================
# Step 7: Packing list to save in loop
## ====================================================================
ls[[i]] <- list("RegularisedCoefs"=reg_coef_tbl,
"bestA"=best_alph,
"bestL"=best_lamb,
"ConfusionMatrx"=reg_cfm,
"AUROC"=reg_auc_sum)

41
1 PA Decline/Til DDV/standardise.R Normal file
View file

@ -0,0 +1,41 @@
## ItMLiHSmar2022
## standardise.R, child script
## Data standardisation, returns list
## Andreas Gammelgaard Damsbo, agdamsbo@clin.au.dk
standardise<-function(train,test,type){
# From:
# https://datascience.stackexchange.com/questions/13971/standardization-normalization-test-data-in-r
sel<-sapply(Xtrain,is.numeric) # Deciding which to stadardise (only numeric)
cnm<-colnames(Xtrain) # Saving column names for ordering
# Subsetting
## Data to treat
train.tr<-train[,sel]
test.tr<-test[,sel]
## Data to save
train.sv<-train[,!sel]
test.sv<-test[,!sel]
# Calculate mean and SD of train data
trainMean <- sapply(train.tr,mean)
trainSd <- sapply(train.tr,sd)
if (type=="c"){
## centered
norm.trainData<-sweep(train.tr, 2L, trainMean) # using the default "-" to subtract mean column-wise
norm.testData<-sweep(test.tr, 2L, trainMean) # using the default "-" to subtract mean column-wise
}
if (type=="cs"){
## centered AND scaled (Z-score standardisation)
norm.trainData<-sweep(sweep(train.tr, 2L, trainMean), 2, trainSd, "/")
norm.testData<-sweep(sweep(test.tr, 2L, trainMean), 2, trainSd, "/")
}
return(list(XtrainSt=cbind(norm.trainData,train.sv)[,cnm], # Reordering columns to original
XtestSt=cbind(norm.testData,test.sv)[,cnm]))
}