dataPreprocess = function(sim_data, type=NULL, knots=NULL, basis=NULL, nci=NULL,control.points=NULL ){
if(type == "2var"){
# #convert the categorical covariates into dummy variables
rec_obj <- recipe(fatigue_mean_intensity_score_diff_sc ~ ., data = data.sim) # made change here!!!
ind_vars <- rec_obj %>%
step_dummy(all_nominal_predictors())
trained_rec <- prep(ind_vars, training = data.sim)
train_data <- bake(trained_rec, new_data = data.sim)
# #Creating the column for interval 1
# int2 = which(X.mat$int.time_X2 == 1)
# int3 = which(X.mat$int.time_X3 == 1)
# int4 = which(X.mat$int.time_X4 == 1)
# int5 = which(X.mat$int.time_X5 == 1)
# int6 = which(X.mat$int.time_X6 == 1)
# int1 = setdiff(1:275, sort(c(int2, int3, int4, int5, int6)))
int.time_X1 = rep(0, nrow(train_data))
int.time_X1[int1] = 1
#selecting the relevant columns to form the design matrix X
# X.mat <- data.frame(int.time_X1,train_data[, c(34:38, 23,24, 39:50)])
X.mat <- data.frame(int.time_X1,train_data[, c(19:23, 1, 2, 7:18)])
#splitting up the MVPA and ABD into the control and treatment groups
MVPA_C <- ifelse(train_data$Randomization == 0, train_data$mvpa_avg_as_is_1952_orig_diff_sc, 0)
MVPA_I <- ifelse(train_data$Randomization == 1, train_data$mvpa_avg_as_is_1952_orig_diff_sc, 0)
ABD_C <- ifelse(train_data$Randomization ==0, train_data$mean_activity_bout_duration10_diff_sc,0)
ABD_I <- ifelse(train_data$Randomization ==1, train_data$mean_activity_bout_duration10_diff_sc,0)
return(list("X.mat" = X.mat, "MVPA_C" = MVPA_C, "MVPA_I" = MVPA_I, "ABD_C" = ABD_C, "ABD_I" = ABD_I))
}else{
# Create variable for mean activity bout duration
## activity threshold = 10 counts/min
sim_data$mean_activity_bout_duration10 <- 1 / sim_data$astp_as_is_10_orig
# Discard missing response values
sim_data <- sim_data[which(!is.na(sim_data[ , c("fatigue_mean_intensity_score")])),]
# Vector of variables to take difference from baseline
vars <- c("fatigue_mean_intensity_score", "fatigue_mean_interference_score",
"mvpa_avg_as_is_1952_orig", "mean_activity_bout_duration10")
# Variable indicating group assignment
group <- "Randomization"
# Pre-process the data
data <- preprocess_data(sim_data, vars, group)
# Vectors of all the covariates for analysis
## response
y <- "fatigue_mean_intensity_score_diff"
y_sc <- "fatigue_mean_intensity_score_diff_sc"
# y <- "fatigue_mean_interference_score_diff"
# y_sc <- "fatigue_mean_interference_score_diff_sc"
## physical activity vars
pa_vars <- c("mvpa_avg_as_is_1952_orig_diff_sc", "mean_activity_bout_duration10_diff_sc")
pa_vars_orig <- c("mvpa_avg_as_is_1952_orig_diff", "mean_activity_bout_duration10_diff")
plot.title <- c("MVPA", "Mean ABD") # titles for plotting
## non-PA covariates
covars <- c("age100", "BMI_obese", "ht_sample", "breast_sample")
covars <- c(paste0(covars, ".C"), paste0(covars, ".I"))
#################################################################################
# Change order of data to be blocks for each group and time period. This makes
# plotting curve estimated in JAGS easier. New data will be control group at M3,
# control group at M6, ..., treatment group at month 6, treatment group at
# month 12.
################################################################################
## Separate data by group (control vs intervention)
indC <- which(data$Randomization == 0)
indI <- which(data$Randomization == 1)
dataC <- data[indC, ]
dataI <- data[indI, ]
### Get indices for time periods for each group
# Month 3
t1.C <- which(data$svy_interval[indC] == 1)
t1.I <- which(data$svy_interval[indI] == 1)
# Month 6
t2.C <- which(data$svy_interval[indC] == 2)
t2.I <- which(data$svy_interval[indI] == 2)
# Month 12
t3.C <- which(data$svy_interval[indC] == 3)
t3.I <- which(data$svy_interval[indI] == 3)
### Reorder response
yC <- as.matrix(dataC[ , y_sc], ncol = 1)
yI <- as.matrix(dataI[ , y_sc], ncol = 1)
colnames(yC) <- colnames(yI) <- y_sc
Y <- as.matrix(c(yC[c(t1.C, t2.C, t3.C)], yI[c(t1.I, t2.I, t3.I)]), ncol = 1)
### Reorder matrix for PA variables
xC <- as.matrix(dataC[ , pa_vars], ncol = length(pa_vars))
xI <- as.matrix(dataI[ , pa_vars], ncol = length(pa_vars))
X <- rbind(xC[c(t1.C, t2.C, t3.C),], xI[c(t1.I, t2.I, t3.I),])
### Reorder matrix for other covariates
zC <- as.matrix(dataC[ , covars[grepl(".C", covars, fixed = TRUE)]]) # made change here!!!!
zI <- as.matrix(dataI[ , covars[grepl(".I", covars, fixed = TRUE)]])
Z <- cbind(rbind(zC[c(t1.C, t2.C, t3.C),], matrix(0, nrow = nrow(zI), ncol = ncol(zC))),
rbind(matrix(0, nrow = nrow(zC), ncol = ncol(zI)), zI[c(t1.I, t2.I, t3.I),]))
## Build a dataset containing variables of interest
int.and.group <- rbind(dataC[c(t1.C, t2.C, t3.C), c("int.time", group)],
dataI[c(t1.I, t2.I, t3.I), c("int.time", group)])
data.sim <- data.frame(Y, Z, int.and.group, X)
colnames(data.sim) <- c(y_sc, colnames(zC), colnames(zI), "int.time", group, pa_vars)
# Change age variables to numeric
# (turns into factor when combined into matrix with other factor variables)
data.sim[, covars[grepl("age", covars)]] <- apply(data.sim[, covars[grepl("age", covars)]], 2, as.numeric)
#convert the categorical covariates into dummy variables
# Formula for model for mean ABD and MVPA on fatigue
k <- knots # number of knots for spline
form <- paste0(y_sc, " ~ -1 + int.time")
for(pa in pa_vars){
form <- paste0(form, " + s(", pa, ", bs = basis, k = ", k, ", by = as.factor(Randomization))")
}
for(param in c(colnames(zC), colnames(zI))){
form <- paste0(form, " + ", param)
}
# Fit model for mean ABD and MVPA on fatigue intensity
fit <- gam(as.formula(form), data = data.sim,
method = 'REML', select = TRUE)
jags.file <- paste(function_directory, "/code/","/test.jags",sep="")
fit.jagam <- jagam(as.formula(form),data=data.sim, file=jags.file,
sp.prior="gamma",diagonalize=TRUE)
# Derive linear predictor matrix
X.mat.full <- predict.gam(fit, type = "lpmatrix")
X.mat <- X.mat.full[ ,1:nci] #covariates
W1_C.mat <- X.mat.full[ ,(nci+1):(nci+control.points-1)] # MVPA for control group for different knot points
W1_I.mat <- X.mat.full[ ,(nci+control.points):(nci+2*control.points-2)] # MVPA for treatment group for different knot points
W2_C.mat <- X.mat.full[ ,(nci+2*control.points-1):(nci+3*control.points-3)] # ABD for control group for different knot points
W2_I.mat <- X.mat.full[ ,(nci+3*control.points-2):(nci+4*control.points-4)] # ABD for treatment group for different knot points
S1 <- cbind(fit$smooth[[1]]$S[[1]], fit$smooth[[1]]$S[[2]]) ## penalty matrices for MVPA control group
S2 <- cbind(fit$smooth[[2]]$S[[1]], fit$smooth[[2]]$S[[2]]) ## penalty matrices for MVPA treatment group
S3 <- cbind(fit$smooth[[3]]$S[[1]], fit$smooth[[3]]$S[[2]]) ## penalty matrices for ABD control group
S4 <- cbind(fit$smooth[[4]]$S[[1]], fit$smooth[[4]]$S[[2]]) ## penalty matrices for ABD treatment group
return(list("X.mat"=X.mat,"W1_C.mat"=W1_C.mat,"W1_I.mat"=W1_I.mat, "W2_C.mat"=W2_C.mat,"W2_I.mat"=W2_I.mat,
"S1"=S1,"S2"=S2,"S3"=S3,"S4"=S4))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.