-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathinitialization.R
62 lines (40 loc) · 2.3 KB
/
initialization.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
##################################################################################################################################
##################################################################################################################################
#load your data here
data<-#read.csv() / read.delim() / read_excel() (require package readxl)
#Name of the columns of the individual indices
patId<-"patId" #labels should be standardized with values from 1 to the number of patient
#Name of the columns of the dose
dose<-"doseStd"
#Name of the columns of the cycle
cycle<-"cycle"
#insert the names of the columns of the grade of the 5 types of toxicity
nams<-c("cutaneous","digestive","genDisorder","hemato","other") #for c("Cutaneous","Digestive","General disorder","Hematologic","Others")
##################################################################################################################################
##################################################################################################################################
library(parallel)
library(rstan)
library(ordinal)
options(mc.cores = 4)
source("src/stanDMfcts.R")
source("src/writeStanCRModel.R")
#Initialize from univariate fixed effect full continuation ratio model, frequentist version for speed
fm1<-clm(as.factor(data[,nams[1]])~1,nominal=~data[,dose]+data[,cycle])
fm2<-clm(as.factor(data[,nams[2]])~1,nominal=~data[,dose]+data[,cycle])
fm3<-clm(as.factor(data[,nams[3]])~1,nominal=~data[,dose]+data[,cycle])
fm4<-clm(as.factor(data[,nams[4]])~1,nominal=~data[,dose]+data[,cycle])
fm5<-clm(as.factor(data[,nams[5]])~1,nominal=~data[,dose]+data[,cycle])
#Write Stan code for univariate model
code<-writeStanCRModel(Y="Y",nY=1,nlvlY=4,X="X",T="T",lvlSlpX=list(1:3),lvlSlpT=list(1:3),sub=TRUE)
outFile<-file("stanCodes/subModel.stan")
writeLines(code$stanFile,outFile,sep="")
close(outFile)
for(i in 1:5){
#get initialization values
initBetaFx<-stdInitBeta(coef(get(paste0("fm",i))))
#Create data for Stan fit
stan_data<-createStanData(data[,nams[i]],data[,dose],data[,cycle],data[,patId],length(initBetaFx))
#fit
assign(paste0("fit",i),stan(file="stanCodes/subModel.stan",data=stan_data,warmup=500,iter=1000,chain=1,init=list(list(betaFx=initBetaFx)),control = list(adapt_delta = 0.99)))
}
save.image("results/saveInits.RData")