/* Joint Model for Y_{1ij} and Y_{2ij}*/

data Dental.joint_glim_est_full;
set Dental.glim_est_full_part Dental.glim_est_full_pard;
run;


proc nlmixed data=Dental.Caldata METHOD=ISAMP noad qpoints=400
tech=NRRIDG seed=6929000;
parms/ data = Dental.joint_glim_est_full;
*bounds 0 < rho_bri12 < 1;
pi = constant(’pi’);


/* binary teeth existence */
uni1 = probnorm(r1);
phi1 = 1.0/sqrt(1 + 3/pi/pi*std_bri1*std_bri1);
B1 = 1.0/phi1*log(sin(pi*uni1*phi1)/sin( phi1*pi*(1-uni1) ) );
xb_bin1= B1 + t_0 + t_age*age + t_sex*sex + t_bmi*bmi+t_smoker*smoker + t_hba1c*hba1c+ t_maxilla*maxilla + t_pov*pov + t_bfl*bfl;
p1 = exp(xb_bin1)/(1 + exp(xb_bin1) );
llik1 = teeth_cal*log( p1 ) + (1 - teeth_cal)*log(1 - p1 );

/* binary disease */
uni2 = probnorm(r2);
phi2 = 1.0/sqrt(1 + 3/pi/pi*std_bri2*std_bri2);
B2 = 1.0/phi2*log(sin(pi*uni2*phi2)/sin( phi2*pi*(1-uni2) ) );
xb_bin2= B2 + d_0 + d_age*age + d_sex*sex + d_bmi*bmi+d_smoker*smoker + d_hba1c*hba1c + d_maxilla*maxilla + d_pov*pov + d_bfl*bfl;
p2 = exp(xb_bin2)/(1 + exp(xb_bin2) );
llik2 = teeth_cal*(d_cal*log(1-p2)+(1 - d_cal)*log(p2));

if d_cal =. then ll = llik1 ;
else ll = llik1 + llik2;

model zzz ~ general(ll);
random r1 r2 ~ normal( [0,0], [1,rho_bri12 , 1] ) subject=id;

ESTIMATE 't_0m' t_0*phi1;
ESTIMATE 't_agem' t_age*phi1;
ESTIMATE 't_sexm' t_sex*phi1;
ESTIMATE 't_bmim' t_bmi*phi1;
ESTIMATE 't_smokerm' t_smoker*phi1;
ESTIMATE 't_hba1cm' t_hba1c*phi1;
ESTIMATE 't_maxillam' t_maxilla*phi1;
ESTIMATE 't_povm' t_pov*phi1;
ESTIMATE 't_bflm' t_bfl*phi1;
ESTIMATE 'd_0m' d_0*phi2;
ESTIMATE 'd_agem' d_age*phi2;
ESTIMATE 'd_sexm' d_sex*phi2;
ESTIMATE 'd_bmim' d_bmi*phi2;
ESTIMATE 'd_smokerm' d_smoker*phi2;
ESTIMATE 'd_hba1cm' d_hba1c*phi2;
ESTIMATE 'd_maxillam' d_maxilla*phi2;
ESTIMATE 'd_povm' d_pov*phi2;
ESTIMATE 'd_bflm' d_bfl*phi2;
estimate 'phi1' phi1;
estimate 'phi2' phi2;
estimate 'rho_bri12' rho_bri12;
run;