R_programs.R
####################################
#Required packages
library(statmod)
library(splines)
library(MASS)
library(nlme)
###################################
mjm.like.aq <- function(par,Dat,knots,BK){
## USAGE
# mjm.like.aq(par,Dat,BK) in conjunction with an
# optimization program, such as optim or nls.
## ARGUMENTS
# par: The values of the following parameters
# beta: a vector of length q (= dim(X)[2]) of regression coefficients to be used in the intercourse model.
# sigma: log standard deviation of the baseline random effect in the intercourse model.
# sigK: log standard deviation of the peakedness random effect in the intercourse model.
# sigY: log standard deviation of the baseline random effect in the pre-ovular model.
# xi: vector of 2 parameters for the effect the intercourse random effects have on the TTP model.
# lambda: vector of max(T) baseline intercepts (i.e., an intercept for each cycle)
# betaT: a vector of length pT (i.e., dim(X_T)[2]) to be used for regression coefficients in the TTP model
# phi: coefficients for the cubic B-splines function.
# rho: a vector coefficients of length pk for the peakedness function.
# cov12: correlation between the baseline and peakedness random effects.
# mu_po,sigma_po: coefficients for the location and scale of the lognormal preovular model.
# Dat: This is a list of all data. Which includes:
# $Y: a vector of intercourse indicators (length L= total number of days observed),
# $ind: a vector of indicators of person number, should be consecutively numbered 1:n (length L),
# $ind_po: a vector of indicators of all cycles being observed (length n=number of subjects),
# $X: intercourse covariate matrix (dimension L x q),
# $X_P: intercourse peakedness covariate matrix (dimension L x pk),
# $day: vector of day relative to ovulation (length L, set to the day of cycle if unknown),
# $day_miss: vector of indicators that the day relative to ovulation known (length L),
# $T: a vector of TTP or Censoring values (length n),
# $CEN: a vector of censoring indicators (length n),
# $X_T: a matrix of covariates for the TTP model (dimension n x pT),
# $nij: a matrix of pre-ovular lengths (dimension n x max(T), set to -1 if unobserved).
# knots: the values of the knots for the cubic b-spline function.
# BK: boundary knots for the cubic b-spline function.
## OUTPUT
# -sum(log(likelihood))
## DETAILS
# This program evaluates a joint model of intercourse and TTP. Pre-ovular length is also
# modeled, since the intercourse model is a function of pre-ovular length and pre-ovular
# length presents informative missingness. The likelihood is evaluated using adaptive
# Gaussian quadrature. First the people with full pre-ovular information are evaluated.
# Then those with missing pre-ovular information are evaluated.
Y <- Dat$Y
ind <- Dat$ind
X <- Dat$X
X_P <- Dat$X_P
X_T <- Dat$X_T
T <- Dat$T
CEN <- Dat$CEN
nij <- Dat$nij
day <- Dat$day
day_miss <- Dat$day_miss
ind_po <- Dat$ind_po
q <- dim(X)[2]
M <- length(T)
tau <- 6
#Defining all parameter values
beta <- par[1:q]
sigma <- exp(par[q+1])+ 1e-10
sigK <- exp(par[q+2])+ 1e-10
sigY <- exp(par[q+3])+ 1e-10
xi <- par[(q+4):(q+5)]
lambda <- exp(par[(q+6):(q+5+tau)])+ 1e-10
betaT <- par[(q+6+tau):(q+5+tau+dim(X_T)[2])]
phi <- par[(q+6+tau+dim(X_T)[2]):(q+5+tau+dim(X_T)[2]+length(knots)+3)]
rho <- par[(q+6+tau+dim(X_T)[2]+length(knots)+3):(q+5+tau+dim(X_T)[2]+length(knots)+3+dim(X_P)[2])]
cov12 <- par[length(par)-2]
cov12 <- pnorm(cov12,0,1)*2-1
if(abs(cov12)==1){cov12 <- sign(cov12)*1 - sign(cov12)*1e-6}
mu_po <- par[length(par)-1]
sig_po <- exp(par[length(par)]) + 1e-10
D <- matrix(c(1,cov12,0,cov12,1,0,0,0,1),3,3)
#Calculating regression values
eta <- X%*%beta
eta_T <- X_T%*%betaT
eta_K <- X_P%*%rho
#First AGQ will be implemented for those will no missing pre-ovular lengths
unq.id <- unique(ind)
no_miss_ind <- unq.id[ind_po==1]
count <- 1
md.est <- NULL
RE.v <- gauss.quad.prob(7,"normal",sigma=1)
t.quad.mat <- trip.quad.prob(RE.v,RE.v,RE.v,7)
re.sum2 <- apply(t.quad.mat$node.mat^2,1,sum)
Nint.mat <- NULL
stop.ind <- 0
if(any(!is.finite(c(beta,sigma,sigK,sigY,xi,lambda,betaT,phi,rho,cov12,mu_po,sig_po)))){stop.ind=1}
for(id in no_miss_ind){
if(count > 1){
if(any(is.na(Nint.mat))){
stop.ind = 1
}
if(any(is.nan(Nint.mat))){
stop.ind = 1
}
if(all(!is.nan(Nint.mat)) & all(!is.na(Nint.mat))){
if(any(Nint.mat==0)){
stop.ind = 1
}
}
}
Nall.mat = 0
if(stop.ind==0){
#Extracting necessary quantities for persion 'id'
tY <- Y[ind==id]
tdelt <- eta[ind==id]
teta_K <- eta_K[ind==id]
tX_P.d <- X_P[ind==id,]
tday <- day[ind==id]-1
teta_T <- eta_T[id]
tTTP <- T[id]
tCEN <- CEN[id]
men.1 <- nij[id,nij[id,]>(-1)]
#Days outside of the boudary knots set BK
tday[tday>BK[2]] <- BK[2]
tday[tday<BK[1]] <- BK[1]
#Calculating the b-splines contribution for person 'id'
basis <- bs(tday,knots=knots,Boundary.knots = BK)
X_B <- cbind(basis[(1:dim(basis)[1]),1:dim(basis)[2]])
teta_B <- X_B%*%phi
#Implementing AGQ. First find the mode of the posterior.
md.opt<-try(optim(c(0,0,0),like.int.gammaj.sh3,tdelt=tdelt,teta_T=teta_T,teta_B=teta_B,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1,method="BFGS",hessian=TRUE),silent=TRUE)
if(class(md.opt)!="try-error"){
hess2 <- md.opt$hessian
bhat2 <- md.opt$par
test <- try(var <- ginv(hess2),silent=TRUE)
if(class(test)=="try-error" | !is.numeric(var) | any(!is.finite(var))){var <- -ginv(diag(3))}
if(!is.numeric(eigen(var)$values)){var <- -ginv(diag(3))}
###If the variance matrix is not proper. Use a different method to calculate##
if(any(eigen(var)$values<0)){
t1 <- fdHess(bhat2,like.int.gammaj.sh3,tdelt=tdelt,teta_T=teta_T,teta_B=teta_B,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1)
hess2 <- t1$Hessian
test <- try(var2 <- ginv(hess2),silent=TRUE)
if(class(test)=="try-error"){var2 <- ginv(-diag(3))}
if(any(eigen(var2)$values<0)){
hess2 <- diag(c(1,1,1))
bhat2 <- c(0,0,0)
}
}
}
#Testing that the AGQ was succesful. If it failed we switch to standard GQ
if(class(md.opt)=="try-error"){
hess2 <- diag(3)
bhat2 <- c(0,0,0)
}
if(any(eigen(ginv(hess2))$values<0) | det(hess2)<0 ){
hess2 <-diag(3)
bhat2 <- c(0,0,0)
}
cho.INV <- try(chol(ginv(hess2)),silent=TRUE)
if(class(cho.INV)=="try-error"){
cho.INV <- diag(3)
bhat2 <- c(0,0,0)
hess2 <- diag(3)
}
#Calculating the shifted and scaled quadrature nodes
btt1 <- t.quad.mat$node.mat %*% cho.INV
btilde <- t(bhat2 + t(btt1))
#Calculating the likelihood contribution
i.like.vals <- apply(btilde,1,function(x){like.int.gammajI.sh3(x,tdelt=tdelt,teta_T=teta_T,teta_B=teta_B,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1)})
i.like.vals[is.nan(i.like.vals)] <- 0
i.like.vals[is.na(i.like.vals)] <- 0
i.like.vals[i.like.vals==Inf] <- 1e100
Nall.mat <- (2*pi)^(3/2)*sum(exp(re.sum2/2)*dmnorm(btilde,rep(0,dim(D)[1]),D)*t.quad.mat$weight.mat*i.like.vals)/sqrt(det(hess2))
if(!is.nan(Nall.mat)){if(!is.na(Nall.mat)){if(Nall.mat==0){Nall.mat <- -1000}}}
}
Nint.mat <- c(Nint.mat,Nall.mat)
count <- count +1
}
#Second AGQ will be implemented for those with a missing pre-ovular length
unq.id <- unique(ind)
miss_ind <- unq.id[ind_po==0]
md.est <- NULL
RE.v <- gauss.quad.prob(7,"normal",sigma=1)
t.quad.mat <- quad.quad.prob(RE.v,RE.v,RE.v,RE.v,7)
re.sum2 <- apply(t.quad.mat$node.mat^2,1,sum)
Nint.mat2 <- NULL
count <- 1
for(id in miss_ind){
if(count > 1){
if(any(is.na(Nint.mat2))){
stop.ind = 1
}
if(any(is.nan(Nint.mat2))){
stop.ind = 1
}
if(all(!is.nan(Nint.mat2)) & all(!is.na(Nint.mat2))){
if(any(Nint.mat2==0)){
stop.ind = 1
}
}
}
Nall.mat = 0
if(stop.ind==0){
#Extracting necessary quantities for persion 'id'
tY <- Y[ind==id]
tdelt <- eta[ind==id]
teta_K <- eta_K[ind==id]
tX_P.d <- X_P[ind==id,]
tday <- day[ind==id]-1
tday_mis<- day_miss[ind==id]
tday[tday_mis==0] <- tday[tday_mis==0]+1
teta_T <- eta_T[id]
tTTP <- T[id]
tCEN <- CEN[id]
men.1 <- nij[id,1:tTTP]
#Implementing AGQ. First find the mode of the posterior.
md.opt<-try(optim(c(0,0,0,mu_po),like.int.gammaj.sh4,tdelt=tdelt,teta_T=teta_T,tday=tday,tday_mis=tday_mis,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,phi=phi,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1,knots=knots,BK=BK,method="BFGS",hessian=TRUE),silent=TRUE)
if(class(md.opt)!="try-error"){
hess2 <- md.opt$hessian
bhat2 <- md.opt$par
test <- try(var <- ginv(hess2),silent=TRUE)
if(class(test)=="try-error" | !is.numeric(var) | any(!is.finite(var))){var <- ginv(-diag(4))}
if(!is.numeric(eigen(var)$values)){var <- -ginv(diag(3))}
###If the variance matrix is not proper. Use a different method to calculate##
if(any(eigen(var)$values<0)){
t1 <- fdHess(bhat2,like.int.gammaj.sh4,tdelt=tdelt,teta_T=teta_T,tday=tday,tday_mis=tday_mis,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,phi=phi,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1,knots=knots,BK=BK)
hess2 <- t1$Hessian
test <- try(var2 <- ginv(hess2),silent=TRUE)
if(class(test)=="try-error"){var2 <- ginv(-diag(4))}
if(any(eigen(var2)$values<0)){
hess2 <- diag(4)
bhat2 <- c(0,0,0,0)
}
}
}
#Testing that the AGQ was succesful. If it failed we switch to standard GQ
if(class(md.opt)=="try-error"){
hess2 <- diag(4)
bhat2 <- c(0,0,0,0)
}
if(any(eigen(ginv(hess2))$values<0)| det(hess2)<0 ){
hess2 <- diag(4)
bhat2 <- c(0,0,0,0)
}
cho.INV <- try(chol(ginv(hess2)),silent=TRUE)
if(class(cho.INV)=="try-error"){
cho.INV <- diag(4)
bhat2 <- c(0,0,0,0)
hess2 <- diag(4)
}
#Calculating the shifted and scaled quadrature nodes
btt1 <- t.quad.mat$node.mat %*% cho.INV
btilde <- t(bhat2 + t(btt1))
btilde2<- cbind(btilde[,1:3],btilde[,4]-(mu_po+sigY*btilde[,3])*exp(sig_po^2/2))
if(any(abs(btilde2)==Inf)){btilde2[abs(btilde2)==Inf] <- sign(btilde2[abs(btilde2)==Inf])*1e16}
#Covariance matrix of btilde2
D2 <- diag(4)
D2[c(1:3),c(1:3)] <- D
D2[4,4] <- sig_po
#Calculating the likelihood contribution
i.like.vals <- apply(btilde,1,function(x){like.int.gammajI.sh4(x,tdelt=tdelt,teta_T=teta_T,tday=tday,tday_mis=tday_mis,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,phi=phi,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1,knots=knots,BK=BK)})
i.like.vals[is.nan(i.like.vals)] <- 0
i.like.vals[is.na(i.like.vals)] <- 0
i.like.vals[i.like.vals==Inf] <- 1e100
Nall.mat <- (2*pi)^(4/2)*sum(exp(re.sum2/2)*dmnorm(btilde2,c(0,0,0,0),D2)*t.quad.mat$weight.mat*i.like.vals)/sqrt(det(hess2))
if(!is.nan(Nall.mat)){if(!is.na(Nall.mat)){if(Nall.mat==0){Nall.mat <- -1000}}}
}
Nint.mat2 <- c(Nint.mat2,Nall.mat)
count <- count +1
}
All_like <- c(Nint.mat,Nint.mat2)
All_like[is.na(All_like)] <- -1000
All_like[is.nan(All_like)] <- -1000
L.All_Like <- All_like
L.All_Like[All_like<=0] <- -1000
L.All_Like[All_like>0] <- log(All_like[All_like>0])
s1 = sum(-L.All_Like)
return(s1)
}
mjm.EB <- function(par,Dat,nodes,knots,BK){
## USAGE
# mjm.EB(par,Dat,nodes,knots,BK): used to get empirical Bayes estimates of the random effects.
## ARGUMENTS
# par: The values of the following parameters
# beta: a vector of length q (= dim(X)[2]) of regression coefficients to be used in the intercourse model.
# sigma: log standard deviation of the baseline random effect in the intercourse model.
# sigK: log standard deviation of the peakedness random effect in the intercourse model.
# sigY: log standard deviation of the baseline random effect in the pre-ovular model.
# xi: vector of 2 parameters for the effect the intercourse random effects have on the TTP model.
# lambda: vector of max(T) baseline intercepts (i.e., an intercept for each cycle)
# betaT: a vector of length pT (i.e., dim(X_T)[2]) to be used for regression coefficients in the TTP model
# phi: coefficients for the cubic B-splines function.
# rho: a vector coefficients of length pk for the peakedness function.
# cov12: correlation between the baseline and peakedness random effects.
# mu_po,sigma_po: coefficients for the location and scale of the lognormal preovular model.
# Dat: This is a list of all data. Which includes:
# $Y: a vector of intercourse indicators (length L= total number of days observed),
# $ind: a vector of indicators of person number, should be consecutively numbered 1:n (length L),
# $ind_po: a vector of indicators of all cycles being observed (length n=number of subjects),
# $X: intercourse covariate matrix (dimension L x q),
# $X_P: intercourse peakedness covariate matrix (dimension L x pk),
# $day: vector of day relative to ovulation (length L, set to the day of cycle if unknown),
# $day_miss: vector of indicators that the day relative to ovulation known (length L),
# $T: a vector of TTP or Censoring values (length n),
# $CEN: a vector of censoring indicators (length n),
# $X_T: a matrix of covariates for the TTP model (dimension n x pT),
# $nij: a matrix of pre-ovular lengths (dimension n x max(T), set to -1 if unobserved).
# nodes: number of quadrature nodes.
# knots: the values of the knots for the cubic b-spline function.
# BK: boundary knots for the cubic b-spline function.
## OUTPUT
# A n x 3 matrix of empirical Bayes Random Effects.
## DETAILS
# This program produces empirical Bayes estimates of the Random Effects.
Y <- Dat$Y
ind <- Dat$ind
X <- Dat$X
X_P <- Dat$X_P
X_T <- Dat$X_T
T <- Dat$T
CEN <- Dat$CEN
nij <- Dat$nij
day <- Dat$day
day_miss <- Dat$day_miss
ind_po <- Dat$ind_po
q <- dim(X)[2]
M <- length(T)
tau <- 6
#Defining all parameter values
beta <- par[1:q]
sigma <- exp(par[q+1])+ 1e-10
sigK <- exp(par[q+2])+ 1e-10
sigY <- exp(par[q+3])+ 1e-10
xi <- par[(q+4):(q+5)]
lambda <- exp(par[(q+6):(q+5+tau)])+ 1e-10
betaT <- par[(q+6+tau):(q+5+tau+dim(X_T)[2])]
phi <- par[(q+6+tau+dim(X_T)[2]):(q+5+tau+dim(X_T)[2]+length(knots)+3)]
rho <- par[(q+6+tau+dim(X_T)[2]+length(knots)+3):(q+5+tau+dim(X_T)[2]+length(knots)+3+dim(X_P)[2])]
cov12 <- par[length(par)-2]
mu_po <- par[length(par)-1]
sig_po <- exp(par[length(par)]) + 1e-10
D <- matrix(c(1,cov12,0,cov12,1,0,0,0,1),3,3)
#Calculating regression values
eta <- X%*%beta
eta_T <- X_T%*%betaT
eta_K <- X_P%*%rho
#First EB will be implemented for those will no missing pre-ovular lengths
unq.id <- unique(ind)
no_miss_ind <- unq.id[ind_po==1]
EB.mat <- NULL
Hess.mat <- NULL
stop.ind <- 0
for(id in no_miss_ind){
#Extracting necessary quantities for persion 'id'
tY <- Y[ind==id]
tdelt <- eta[ind==id]
teta_K <- eta_K[ind==id]
tX_P.d <- X_P[ind==id,]
tday <- day[ind==id]-1
teta_T <- eta_T[id]
tTTP <- T[id]
tCEN <- CEN[id]
men.1 <- nij[id,nij[id,]>(-1)]
#Days outside of the boudary knots set BK
tday[tday>BK[2]] <- BK[2]
tday[tday<BK[1]] <- BK[1]
#Calculating the b-splines contribution for person 'id'
basis <- bs(tday,knots=knots,Boundary.knots = BK)
X_B <- cbind(basis[(1:dim(basis)[1]),1:dim(basis)[2]])
teta_B <- X_B%*%phi
#Implementing EB.
md.opt<-optim(c(0,0,0),like.int.gammaj.sh3,tdelt=tdelt,teta_T=teta_T,teta_B=teta_B,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1,method="BFGS",hessian=TRUE)
hess2 <- md.opt$hessian
bhat2 <- md.opt$par
var <- ginv(hess2)
EB.mat <- rbind(EB.mat, bhat2)
Hess.mat <- rbind(Hess.mat,c(var))
}
#Second EB will be implemented for those with a missing pre-ovular length
unq.id <- unique(ind)
miss_ind <- unq.id[ind_po==0]
EB.mat2 <- NULL
Hess.mat2<- NULL
for(id in miss_ind){
#Extracting necessary quantities for persion 'id'
tY <- Y[ind==id]
tdelt <- eta[ind==id]
teta_K <- eta_K[ind==id]
tX_P.d <- X_P[ind==id,]
tday <- day[ind==id]-1
tday_mis<- day_miss[ind==id]
teta_T <- eta_T[id]
tTTP <- T[id]
tCEN <- CEN[id]
men.1 <- nij[id,1:tTTP]
#Implementing EB.
md.opt<- optim(c(0,0,0,mu_po),like.int.gammaj.sh4,tdelt=tdelt,teta_T=teta_T,tday=tday,tday_mis=tday_mis,teta_K = teta_K,tTTP=tTTP,tCEN=tCEN,tY=tY,D=D,mu_po=mu_po,phi=phi,sig_po=sig_po,lambda=lambda,xi=xi,sigma=sigma,sigK=sigK,sigY=sigY,men.1=men.1,knots=knots,BK=BK,method="BFGS",hessian=TRUE)
hess2 <- md.opt$hessian
bhat2 <- md.opt$par
var <- ginv(hess2)
EB.mat2 <- rbind(EB.mat2, bhat2[1:3])
Hess.mat2 <- rbind(Hess.mat,c(var[c(1:3),c(1:3)]))
}
Fin_EB.mat <- NULL
Fin_Var.mat <- NULL
count_miss <- 1
count <- 1
#Re-ordering the EB estimates.
for(j in unq.id){
if(ind_po[unq.id==j]==1){
Fin_EB.mat <- rbind(Fin_EB.mat,EB.mat[count,])
Fin_Var.mat<- rbind(Fin_Var.mat,Hess.mat[count,])
count <- count +1
}
if(ind_po[unq.id==j]==0){
Fin_EB.mat <- rbind(Fin_EB.mat,EB.mat2[count_miss,])
Fin_Var.mat<- rbind(Fin_Var.mat,Hess.mat2[count_miss,])
count_miss <- count_miss +1
}
}
Res <- list(Fin_EB.mat=Fin_EB.mat,Fin_Var.mat=Fin_Var.mat)
return(Res)
}
mjm.pred.probs <- function(par,Dat,knots,BK){
## USAGE
# mjm.pred.probs(par,Dat,nodes,knots,BK): used to get estimates of the predicted intercourse and survival probabilities.
Y <- Dat$Y
ind <- Dat$ind
X <- Dat$X
X_P <- Dat$X_P
X_T <- Dat$X_T
re_mat <- Dat$re_mat
unq.id <- unique(ind)
q <- dim(X)[2]
M <- length(T)
tau <- 6
#Defining all parameter values
beta <- par[1:q]
sigma <- exp(par[q+1])+ 1e-10
sigK <- exp(par[q+2])+ 1e-10
sigY <- exp(par[q+3])+ 1e-10
xi <- par[(q+4):(q+5)]
lambda <- exp(par[(q+6):(q+5+tau)])+ 1e-10
betaT <- par[(q+6+tau):(q+5+tau+dim(X_T)[2])]
phi <- par[(q+6+tau+dim(X_T)[2]):(q+5+tau+dim(X_T)[2]+length(knots)+3)]
rho <- par[(q+6+tau+dim(X_T)[2]+length(knots)+3):(q+5+tau+dim(X_T)[2]+length(knots)+3+dim(X_P)[2])]
cov12 <- par[length(par)-2]
mu_po <- par[length(par)-1]
sig_po <- exp(par[length(par)])
D <- matrix(c(1,cov12,0,cov12,1,0,0,0,1),3,3)
#Calculating regression values
eta_T <- X_T%*%betaT
eta_K <- X_P%*%rho
#Estimated Day
day <- rep(0,length(X[,1]))
for(j in unq.id){
t_day <- day[ind==j]
t_cyc <- X_P[ind==j,1]
t_re <- re_mat[j,]
mean_po<- exp(log(max(c(mu_po+sigY*t_re[3],0.001)))+sig_po^2/2)
for(k in unique(t_cyc)){
tt_day <- t_day[t_cyc==k]
tt_day <- seq(1,length(tt_day)) - mean_po
t_day[t_cyc==k] <- tt_day
}
day[ind==j] <- t_day
}
int.prob.mat <- NULL
surv.prob.mat <- NULL
for(id in unq.id){
#Extracting necessary quantities for persion 'id'
tY <- Y[ind==id]
tX <- X[ind==id,]
teta_K <- eta_K[ind==id]
tday <- day[ind==id]
teta_T <- eta_T[id]
t_re <- re_mat[id,]
#Days outside of the boudary knots set BK
tday[tday>BK[2]] <- BK[2]
tday[tday<BK[1]] <- BK[1]
#Calculating the b-splines contribution for person 'id'
basis <- bs(tday,knots=knots,Boundary.knots = BK)
X_B <- cbind(basis[(1:dim(basis)[1]),1:dim(basis)[2]])
teta_B <- X_B%*%phi
#Implementing EB.
est_probs <- int_surv_prob(t_re,tdelt,teta_T,teta_B,teta_K,tX,beta,lambda,xi,sigma,sigK,tday)
int.prob.mat <- rbind(int.prob.mat,cbind(rep(id,length(tX[,1])),est_probs$int_prob))
surv.prob.mat<- rbind(surv.prob.mat,c(id,est_probs$cond_surv))
}
Res_list <- list(int.prob.mat= int.prob.mat, surv.prob.mat= surv.prob.mat)
return(Res_list)
}
########################################################################################
trip.quad.prob<- function(RE,RE2,RE3,nodes=20){
node.mat=NULL
weight.mat=NULL
for(i in 1:nodes){
for(j in 1:nodes){
node.mat = rbind(node.mat,cbind(rep(RE$nodes[i],nodes),rep(RE2$nodes[j],nodes),RE3$nodes))
weight.mat = c(weight.mat,rep(RE$weights[i],nodes)*rep(RE2$weights[j],nodes)*RE3$weights)
}}
return(list(node.mat=node.mat,weight.mat=weight.mat))
}
quad.quad.prob<- function(RE,RE2,RE3,RE4,nodes=20){
node.mat=NULL
weight.mat=NULL
for(i in 1:nodes){
for(j in 1:nodes){
for(k in 1:nodes){
node.mat = rbind(node.mat,cbind(rep(RE$nodes[i],nodes),rep(RE2$nodes[j],nodes),rep(RE3$nodes[k],nodes),RE4$nodes))
weight.mat = c(weight.mat,rep(RE$weights[i],nodes)*rep(RE2$weights[j],nodes)*rep(RE3$weights[k],nodes)*RE4$weights)
}}}
return(list(node.mat=node.mat,weight.mat=weight.mat))
}
dmnorm <- function (x, mean = rep(0, d), varcov, log = FALSE)
{
d <- if (is.matrix(varcov))
ncol(varcov)
else 1
if (d == 1)
return(dnorm(x, mean, sqrt(varcov), log = log))
x <- if (is.vector(x))
matrix(x, 1, d)
else data.matrix(x)
if (is.vector(mean))
mean <- outer(rep(1, nrow(x)), mean)
if (is.matrix(mean) && (nrow(mean) != nrow(x) || ncol(mean) !=
ncol(x)))
stop("mismatch of dimensions of 'x' and 'mean'")
if (is.vector(mean))
mean <- outer(rep(1, nrow(x)), mean)
X <- t(x - mean)
conc <- pd.solve(varcov, log.det = TRUE)
Q <- apply((conc %*% X) * X, 2, sum)
log.det <- attr(conc, "log.det")
logPDF <- as.vector(Q + d * logb(2 * pi) + log.det)/(-2)
if (log)
logPDF
else exp(logPDF)
}
pd.solve <- function (x, silent = FALSE, log.det = FALSE)
{
if (is.null(x))
return(NULL)
if (any(is.na(x))) {
if (silent)
return(NULL)
else stop("NA's in x")
}
if (length(x) == 1)
x <- as.matrix(x)
if (!is.matrix(x)) {
if (silent)
return(NULL)
else stop("x is not a matrix")
}
if (max(abs(x - t(x))) > .Machine$double.eps) {
if (silent)
return(NULL)
else stop("x appears to be not symmetric")
}
x <- (x + t(x))/2
u <- try(chol(x, pivot = FALSE), silent = silent)
if (class(u) == "try-error") {
if (silent)
return(NULL)
else stop("x appears to be not positive definite")
}
inv <- chol2inv(u)
if (log.det)
attr(inv, "log.det") <- 2 * sum(log(diag(u)))
return(inv)
}
like.int.gammaj.sh3 <-function(re,tdelt,teta_T,teta_B,teta_K,tTTP,tCEN,tY,D,mu_po,sig_po,lambda,xi,sigma,sigK,sigY,men.1){
#Function that is proportional to the posterior distribution for those with no missing PO lengths
RE.mat.F<- sigma*re[1] + teta_B*exp(sigK*re[2]+teta_K)
exp.por <- exp(tdelt + RE.mat.F)
cond.prob<- rep(1,length(exp.por))
cond.prob[is.finite(exp.por)]<- exp.por[is.finite(exp.por)]/(1+exp.por[is.finite(exp.por)])
like.int <- dbinom(tY,1,cond.prob)
like.int[is.nan(like.int)] <- 0
like.int[is.na(like.int)] <- 0
L_int_like <- like.int
L_int_like[like.int==0] <- -1000
L_int_like[like.int>0] <- log(like.int[like.int>0])
like.I <- sum(L_int_like)
l.RE.mat <- lambda*exp(xi[1]*re[1]+xi[2]*re[2]+teta_T)
l.RE.mat <- c(l.RE.mat,l.RE.mat[6])
hazT <- sum(l.RE.mat[1:tTTP])
hazTm1 <- 0
if(tTTP > 1){hazTm1 <- sum(l.RE.mat[1:(tTTP-1)])}
like.sl <- (1-tCEN)*exp(-hazTm1) + (2*tCEN-1)*exp(-hazT)
like.sl[is.nan(like.sl)] <- 0
like.sl[is.na(like.sl)] <- 0
L_SL_like <- like.sl
L_SL_like[like.sl==0] <- -1000
L_SL_like[like.sl>0] <- log(like.sl)
like.RE <- dmnorm(re,rep(0,length(re)),D)
like.RE[is.nan(like.RE)] <- 0
like.RE[is.na(like.RE)] <- 0
L_RE_like <- like.RE
L_RE_like[like.RE==0] <- -1000
L_RE_like[like.RE>0] <- log(like.RE)
like.Y <- dlnorm(men.1[men.1>0],log(min(max(c(mu_po+sigY*re[3],0.001)),1e10)),sig_po)
like.Y[is.nan(like.Y)] <- 0
like.Y[is.na(like.Y)] <- 0
L_MCL_like <- like.Y
L_MCL_like[like.Y ==0] <- -1000
L_MCL_like[like.Y>0] <- log(like.Y[like.Y>0])
like.IY <- sum(L_MCL_like)
llike <- -(like.I + L_SL_like + L_RE_like + like.IY)
llike
}
like.int.gammajI.sh3 <-function(re,tdelt,teta_T,teta_B,teta_K,tTTP,tCEN,tY,D,mu_po,sig_po,lambda,xi,sigma,sigK,sigY,men.1){
#Evaluating the likelihood conditional on the unknown values
RE.mat.F<- sigma*re[1] + teta_B*exp(sigK*re[2]+teta_K)
exp.por <- exp(tdelt + RE.mat.F)
cond.prob<- rep(1,length(exp.por))
cond.prob[is.finite(exp.por)]<- exp.por[is.finite(exp.por)]/(1+exp.por[is.finite(exp.por)])
like.int <- dbinom(tY,1,cond.prob)
like.int[is.nan(like.int)] <- 0
like.int[is.na(like.int)] <- 0
L_int_like <- like.int
L_int_like[like.int==0] <- -1000
L_int_like[like.int>0] <- log(like.int[like.int>0])
like.I <- sum(L_int_like)
l.RE.mat <- lambda*exp(xi[1]*re[1]+xi[2]*re[2]+teta_T)
l.RE.mat <- c(l.RE.mat,l.RE.mat[6])
hazT <- sum(l.RE.mat[1:tTTP])
hazTm1 <- 0
if(tTTP > 1){hazTm1 <- sum(l.RE.mat[1:(tTTP-1)])}
like.sl <- (1-tCEN)*exp(-hazTm1) + (2*tCEN-1)*exp(-hazT)
like.sl[is.nan(like.sl)] <- 0
like.sl[is.na(like.sl)] <- 0
L_SL_like <- like.sl
L_SL_like[like.sl==0] <- -1000
L_SL_like[like.sl>0] <- log(like.sl)
like.Y <- dlnorm(men.1[men.1>0],log(min(max(c(mu_po+sigY*re[3],0.001)),1e10)),sig_po)
like.Y[is.nan(like.Y)] <- 0
like.Y[is.na(like.Y)] <- 0
L_MCL_like <- like.Y
L_MCL_like[like.Y ==0] <- -1000
L_MCL_like[like.Y>0] <- log(like.Y[like.Y>0])
like.IY <- sum(L_MCL_like)
llike <- -(like.I + L_SL_like + like.IY)
exp(-llike)
}
like.int.gammaj.sh4 <-function(re,tdelt,teta_T,tday,tday_mis,teta_K,tTTP,tCEN,tY,D,mu_po,sig_po,phi,lambda,xi,sigma,sigK,sigY,men.1,knots,BK){
#Function that is proportional to the posterior distribution for those with missing PO lengths
men.1[men.1==0] <- exp(re[4])
tday[tday_mis==0] <- tday[tday_mis==0]-exp(re[4])
#Days outside of the boudary knots set BK
tday[tday>BK[2]] <- BK[2]
tday[tday<BK[1]] <- BK[1]
basis <- bs(tday,knots=knots,Boundary.knots = BK)
X_B <- cbind(basis[(1:dim(basis)[1]),1:dim(basis)[2]])
teta_B <- X_B%*%phi
RE.mat.F<- sigma*re[1] + teta_B*exp(sigK*re[2]+teta_K)
exp.por <- exp(tdelt + RE.mat.F)
cond.prob<- rep(1,length(exp.por))
cond.prob[is.finite(exp.por)]<- exp.por[is.finite(exp.por)]/(1+exp.por[is.finite(exp.por)])
like.int <- dbinom(tY,1,cond.prob)
like.int[is.nan(like.int)] <- 0
like.int[is.na(like.int)] <- 0
L_int_like <- like.int
L_int_like[like.int==0] <- -1000
L_int_like[like.int>0] <- log(like.int[like.int>0])
like.I <- sum(L_int_like)
l.RE.mat <- lambda*exp(xi[1]*re[1]+xi[2]*re[2]+teta_T)
l.RE.mat <- c(l.RE.mat,l.RE.mat[6])
hazT <- sum(l.RE.mat[1:tTTP])
hazTm1 <- 0
if(tTTP > 1){hazTm1 <- sum(l.RE.mat[1:(tTTP-1)])}
like.sl <- (1-tCEN)*exp(-hazTm1) + (2*tCEN-1)*exp(-hazT)
like.sl[is.nan(like.sl)] <- 0
like.sl[is.na(like.sl)] <- 0
L_SL_like <- like.sl
L_SL_like[like.sl==0] <- -1000
L_SL_like[like.sl>0] <- log(like.sl)
like.RE <- dmnorm(re[1:3],rep(0,dim(D)[1]),D)
like.RE[is.nan(like.RE)] <- 0
like.RE[is.na(like.RE)] <- 0
L_RE_like <- like.RE
L_RE_like[like.RE==0] <- -1000
L_RE_like[like.RE>0] <- log(like.RE)
like.Y <- dlnorm(men.1[men.1>0],log(min(max(c(mu_po+sigY*re[3],0.001)),1e10)),sig_po)
like.Y[is.nan(like.Y)] <- 0
like.Y[is.na(like.Y)] <- 0
L_MCL_like <- like.Y
L_MCL_like[like.Y ==0] <- -1000
L_MCL_like[like.Y>0] <- log(like.Y[like.Y>0])
like.IY <- sum(L_MCL_like)
llike <- -(like.I + L_SL_like + L_RE_like + like.IY)
llike
}
like.int.gammajI.sh4 <-function(re,tdelt,teta_T,tday,tday_mis,teta_K,tTTP,tCEN,tY,D,mu_po,sig_po,phi,lambda,xi,sigma,sigK,sigY,men.1,knots,BK){
#Evaluating the likelihood conditional on the unknown values (RE and unknown PO length)
men.1[men.1==0] <- exp(re[4])
tday[tday_mis==0] <- tday[tday_mis==0]-exp(re[4])
tday[tday>BK[2]] <- BK[2]
tday[tday<BK[1]] <- BK[1]
basis <- bs(tday,knots=knots,Boundary.knots = BK)
X_B <- cbind(basis[(1:dim(basis)[1]),1:dim(basis)[2]])
teta_B <- X_B%*%phi
RE.mat.F<- sigma*re[1] + teta_B*exp(sigK*re[2]+teta_K)
exp.por <- exp(tdelt + RE.mat.F)
cond.prob<- rep(1,length(exp.por))
cond.prob[is.finite(exp.por)]<- exp.por[is.finite(exp.por)]/(1+exp.por[is.finite(exp.por)])
like.int <- dbinom(tY,1,cond.prob)
like.int[is.nan(like.int)] <- 0
like.int[is.na(like.int)] <- 0
like.I <- sum(log(like.int))
L_int_like <- like.int
L_int_like[like.int==0] <- -1000
L_int_like[like.int>0] <- log(like.int[like.int>0])
like.I <- sum(L_int_like)
l.RE.mat <- lambda*exp(xi[1]*re[1]+xi[2]*re[2]+teta_T)
l.RE.mat <- c(l.RE.mat,l.RE.mat[6])
hazT <- sum(l.RE.mat[1:tTTP])
hazTm1 <- 0
if(tTTP > 1){hazTm1 <- sum(l.RE.mat[1:(tTTP-1)])}
like.sl <- (1-tCEN)*exp(-hazTm1) + (2*tCEN-1)*exp(-hazT)
like.sl[is.nan(like.sl)] <- 0
like.sl[is.na(like.sl)] <- 0
L_SL_like<- like.sl
L_SL_like[like.sl==0] <- -1000
L_SL_like[like.sl>0] <- log(like.sl)
like.Y <- dlnorm(men.1[men.1>0],log(min(max(c(mu_po+sigY*re[3],0.001)),1e10)),sig_po)
like.Y[is.nan(like.Y)] <- 0
like.Y[is.na(like.Y)] <- 0
L_MCL_like <- like.Y
L_MCL_like[like.Y ==0] <- -1000
L_MCL_like[like.Y>0] <- log(like.Y[like.Y>0])
like.IY <- sum(L_MCL_like)
llike <- -(like.I + L_SL_like + like.IY)
exp(-llike)
}
int_surv_prob <-function(re,tdelt,teta_T,teta_B,teta_K,tX,beta,lambda,xi,sigma,sigK,tday){
#Calculating intercourse probabilities taking into account the unknown lag
q <- length(beta)
RE.mat.F <- sigma*re[1] + teta_B*exp(sigK*re[2]+teta_K)
leng_prob <- length(tX[,1])
lag <- 0
int_prob <- NULL
for(si in 1:leng_prob){
tdelt_0 <- sum(beta[1:c(q-1)]*tX[si,1:c(q-1)])
exp.por_0 <- exp(tdelt_0 + RE.mat.F[si,])/(1+exp(tdelt_0 + RE.mat.F[si,]))
exp.por_1 <- exp(tdelt_0 + beta[q] + RE.mat.F[si,])/(1+exp(tdelt_0 + beta[q] + RE.mat.F[si,]))
test_prob <- lag*exp.por_1 + (1-lag)*exp.por_0
lag <- test_prob
int_prob <- c(int_prob,test_prob)
}
#estimating the conditional survival function
FW_probs <- int_prob[tday<=1 & tday>=(-6)]
if(length(int_prob[tday<=1 & tday>=(-6)])==0){
FW_probs <- max(int_prob)
}
z.int.prob <- (1-mean(FW_probs))^8
l.RE.mat <- (1-z.int.prob)*lambda*exp(xi[1]*re[1]+xi[2]*re[2]+teta_T)
cum_hazT <- cumsum(l.RE.mat)
surv <- exp(-cum_hazT)
cond_surv<- surv/surv[1]
est_probs <- list(int_prob=int_prob,cond_surv=cond_surv)
return(est_probs)
}
BACK TO TOP