model-for-jags-NIOP_BX-IOP_SURG.txt

model {





###PRIORS FOR LATENT CLASS MODEL

#flat prior on probability of latent class membership for all subjects
#(rho in pdf)
p_eta ~ dbeta(1,1)


###PRIORS FOR MIXED MODEL
#model correlated random effects distribution

for (index in 1:d_Z) {
	xi[index]~dunif(0,100) #scale parameter, same across classes
	for(k in 1:K){
		mu_raw[index,k]~dnorm(0, 0.01) 
		mu[index,k]<-xi[index] *mu_raw[index,k]}  }
		
for(k in 1:K){ #save this iteration of mu with clearer labels. This is not used in model, just to get posteriors more easily.
	mu_int[k] <- mu[1,k]  
	mu_slope[k] <- mu[2,k]} 

#same covariance matrix (Sigma_B_raw) across latent classes
Tau_B_raw ~ dwish(I_d_Z[,], (d_Z+1))  #this is unscaled covariance matrix
Sigma_B_raw[1:d_Z, 1:d_Z] <- inverse(Tau_B_raw[1:d_Z, 1:d_Z])	
for (index in 1:d_Z){
		sigma[index]<-xi[index]*sqrt(Sigma_B_raw[index,index]) } #take into account scaling when saving elements of the covariance matrix

sigma_int <- sigma[1] # again, this is to track elements of the cov matrix with easier labels
sigma_slope <- sigma[2] 

rho_int_slope <- Sigma_B_raw[1,2] / sqrt(Sigma_B_raw[1,1] * Sigma_B_raw[2,2])
cov_int_slope <- rho_int_slope * sigma_int * sigma_slope * xi[1] * xi[2] #again, take into account scaling (xi) when saving the covariance  


##residual variance, independent of correlated random effects, same across classes
sigma_res ~ dunif(0, 1) 
tau_res <- pow(sigma_res,-2)

##fixed effects
for(index in 1:d_X){
	beta[index] ~ dnorm(0,0.01)}


###PRIORS FOR OUTCOME MODEL & BIOPSY DATA
#Add 1 because last element in gamma is coefficients for class membership eta=1
for(index in 1:(d_V_RC+1)){gamma_RC[index] ~ dnorm(0,0.01)}
for(index in 1:(d_W_SURG+2)){omega_SURG[index] ~ dnorm(0,0.01)} #+2 because includes an interaction; IOP_SURG


###LIKELIHOOD

##latent variable for true cancer state
for(i in 1:n_eta_known){
	eta_data[i] ~ dbern(p_eta)
	eta[i] <- eta_data[i] + 1} #this is for those with path reports from SURG, eta known 

for(i in (n_eta_known+1):n){
	eta_hat[(i-n_eta_known)] ~ dbern(p_eta)
	eta[i] <- eta_hat[(i-n_eta_known)] + 1}  #for those without SURG

##linear mixed effects model for PSA 
#generate random intercept and slope for individual given latent class
#in pdf, B_raw is b_i with a 'u' shape on top. (Official term??)
for (i in 1:n) {
	B_raw[i,1:d_Z] ~ dmnorm(mu_raw[1:d_Z,eta[i]], Tau_B_raw[1:d_Z, 1:d_Z])
	for(index in 1:d_Z){b_vec[i,index] <- xi[index]*B_raw[i,index]} }

#fit LME
for(j in 1:n_obs_psa){ 
	mu_obs_psa[j] <- inprod(b_vec[subj_psa[j],1:d_Z], Z[j,1:d_Z])  + (beta[1]*X[j,1:d_X]) 
	Y[j] ~ dnorm(mu_obs_psa[j], tau_res) }


##all biopsy data

#logistic regression for reclassification 	
for(j in 1:n_rc){
	logit(p_rc[j]) <-inprod(gamma_RC[1:d_V_RC], V_RC[j,1:d_V_RC]) + gamma_RC[(d_V_RC+1)]*equals(eta[subj_rc[j]],2) 
	RC[j] ~ dbern(p_rc[j])
}


#logistic regression for surgery (IOP_SURG)
for(j in 1:n_surg){
	logit(p_surg[j]) <-inprod(omega_SURG[1:d_W_SURG], W_SURG[j,1:d_W_SURG]) + omega_SURG[(d_W_SURG+1)]*equals(eta[subj_surg[j]],2)  + omega_SURG[(d_W_SURG+2)]*equals(eta[subj_surg[j]],2)*W_SURG[j,d_W_SURG]  
	SURG[j] ~ dbern(p_surg[j])
}

	
}