***********************************************************************
Purpose: 	Use of PROC NLMIXED to fit likelihood to address both associated measurement error and non-detectable exposures.


Input data: data1 
A simulated dataset (beta=0.2) from a single group of subjects. Each subject has 1 response measurement (resp_ind=1) and between 2 & 6 exposure measurements (resp_ind=0). Approximately 15% of the exposure measurements fall  below a limit of detection (in this example LOD = -1.07 on the natural log scale).

Data dictionary:
subj = Subject ID
n_obs = Record number for a given subject
resp_ind = indicates the record that contains the response
		measurement (ycen)
cens_ind = indicates whether the exposure measurement (resp_ind=0)
is observed (cens_ind=0) or fell below the LOD (cens_ind=1)
ycen = response measurement when resp_ind=1
     = exposure measurement on the natural log scale when resp_ind=0
mnxcen = the mean of the exposure measurements on the original scale
age = covariate age in years
smoke = covariate indicating whether subject is a smoker

Likelihood function is based on a simplified version of equation 5 (drop 'g')
		   with f* from equation 6 replacing f(R_i|delta_i,C_i,Theta) in eq 5.
**********************************************************************;

/*** OBTAIN STARTING VALUES ***;
 *** INSERT IN THE PARMS STATEMENT OF NLMIXED ***;

proc mixed data=out.data1 method=ml asycov;
  where resp_ind=0;
  class subj;
  model ycen = /s;
  random intercept/type=un subject=subj;
run;
proc mixed data=out.data1 method=ml;
  where resp_ind=1;
  model ycen = mnxcen age smoke/s;
run;
*/

proc nlmixed data=out.data1 tech=trureg maxtime=600;
parms sigsqb=1.9 sigsqw=1.6 mu=1 alpha=7 beta=0.1 sigsq= 13 gam1=-0.23 gam2=2.5;
bounds sigsqb sigsqw sigsq >=0;

pi=2*arsin(1);

muy_i = mu + b_i;
mux_i = exp(muy_i + sigsqw/2);
mur_i = alpha + beta*mux_i + gam1*age + gam2*smoke;

loglike =
  (resp_ind)*(-.5*((ycen-mur_i)/sqrt(sigsq))**2 - 
  		log(sqrt(2*pi*sigsq))) +
  (1-resp_ind)*(1-cens_ind)*(-.5*((ycen-muy_i)/sqrt(sigsqw))**2 - 
  		log(sqrt(2*pi*sigsqw))) +
  (1-resp_ind)*(cens_ind)*log((probnorm((ycen - muy_i)/sqrt(sigsqw))));

model ycen ~ general(loglike);
random b_i  ~ normal(0,sigsqb) subject=subj;
ods output parameterestimates=nlmparm convergencestatus=rc;
run;