Oxford.Anal.R
#####################
# Oxford Data Analysis #
# #
####################
library(MASS)
library(nlme)
library(splines)
day_level <- read.csv(file="Oxford_train_day_level.csv")
women_level <- read.csv(file="Oxford_train_women_level.csv")
Y <- day_level[,1] #Intercourse indicators
ind <- day_level[,2] #Subject numbers
X <-as.matrix(day_level[,3:11]) #Day-level covariates for mean
X_P <- as.matrix(day_level[,12:13]) #Day-level covariates for peak
day <- day_level[,14] #Day relative to ovulation (if known)
day_miss <- day_level[,15] #Indicator that the day is known
X_T <- as.matrix(women_level[,1:2]) #TTP covariates
T <- women_level[,3] #TTP
CEN <- women_level[,4] #Censoring indicator
nij <- as.matrix(women_level[,5:11]) #Length of the follicular phase
ind_po <- women_level[,12] #Indicator the length of the follicular is known.
Dat <- list(Y=Y,ind=ind,X=X,X_P=X_P,X_T=X_T,T=T,CEN=CEN,nij=nij,day=day,day_miss=day_miss,ind_po=ind_po)
knots<- c(-13,-10,-7,-4,-2,1,4,7,10,13)
BK <- c(-20,50)
day[day>BK[2]] <- BK[2]
day[day<BK[1]] <- BK[1]
basis <- bs(day,knots=knots,Boundary.knots = BK)
X_B <- cbind(basis[(1:dim(basis)[1]),1:dim(basis)[2]])
beta <- c(-0.93103565, -0.02366856, -0.22017035, 0.17624416, -2.13985499, 0.08124901, -0.13758868, -0.06379638,-0.220662245)
phi <- c( -0.22089327,-0.033253067, 0.13958528, 0.06346478, 0.58637077, 0.43352161, 1.55697622, -0.41935904, -0.37448738,-0.327740735, -1.47752073, 1.33791314, -1.33968068)
rho <- c(0.12325318, -0.02374317)
mu_po <- c(16.008)
sd_po <- log(0.1804)
sigma <- c(-0.19)
sigK <- c(-0.44)
sigY <- c(0.9)
xi <- c(0.2644,0.33877)
lambda<- c(rep(-1.821,6))
betaT <- c(-0.284636,0.6301)
cov12 <- -0.955
par <- c(beta,sigma,sigK,sigY,xi,lambda,betaT,phi,rho,cov12,mu_po,sd_po)
source("R_programs.R")
filename = paste("Oxford_results.csv")
filename
like.res <- optim(par,mjm.like.aq,Dat=Dat,knots=knots,BK=BK,hessian=TRUE)
####Transforming the parameters####
q <- dim(X)[2]
tau <- 6
par_est <- like.res$par
beta <- par_est[1:q]
sig_b1 <- exp(par_est[q+1])
sig_b2 <- exp(par_est[q+2])
sig_b3 <- exp(par_est[q+3])
xi <- par_est[(q+4):(q+5)]
gamma_j <- par_est[(q+6):(q+5+tau)]
gamma <- par_est[(q+6+tau):(q+5+tau+dim(X_T)[2])]
B_splines <- par_est[(q+6+tau+dim(X_T)[2]):(q+5+tau+dim(X_T)[2]+length(knots)+3)]
phi <- par_est[(q+6+tau+dim(X_T)[2]+length(knots)+3):(q+5+tau+dim(X_T)[2]+length(knots)+3+dim(X_P)[2])]
cov_b1b2<- par_est[length(par_est)-2]
cov_b1b2<- pnorm(cov_b1b2,0,1)*2-1
mu_po <- par_est[length(par_est)-1]
sig_po <- exp(par_est[length(par_est)])
var1 <- diag(ginv(like.res$hessian))
####Estimating the variance of the transformed parameters via the Delta Method####
var1[(q+1):(q+3)] <- c(sig_b1,sig_b2,sig_b3)^2*var1[(q+1):(q+3)]
var1[length(par_est)-2] <- (dnorm(cov_b1b2,0,1)*2)*var1[length(par_est)-2]
var1[length(par_est)] <- sig_po^2*var1[length(par_est)]
trans_par_est <- c(beta, sig_b1, sig_b2, sig_b3,xi,gamma_j,gamma, B_splines,phi,cov_b1b2,mu_po,sd_po)
est <- cbind(trans_par_est,var1,rep(like.res$value,length(trans_par_est)))
rownames(est) <- c(colnames(X),"Sig_b1","Sig_b2","Sig_b3","xi_1","xi_2",rep("gamma",length(gamma_j)),colnames(X_T),rep("B_splines",length(B_splines)),"phi_1","phi_2","cov_b1b2","mu_po","sd_po")
colnames(est) <- c("Estimate","Variance","-Log(like)")
est
write.csv(est,file=filename)
BACK TO TOP