title1 "OPEN Study: Regression Calibration using Proc Mixed";
title2 "Nutrient = Log Protein";
title3 "Main Instruments = FFQ & 24HR, Reference Instrument = UN";

*libname lib "h:\documents\kipnis\open\newdata";
libname lib "C:\pctex_Raymond_Carroll_directory\papers\kipnis.protein.directory\open_study_paper";

PROC IMPORT OUT= WORK.fake_open_data 
            DATAFILE= "C:\pctex_Raymond_Carroll_directory\papers\kipnis.
protein.directory\open_study_paper\fake_OPENdata.txt" 
            DBMS=TAB REPLACE;
     GETNAMES=YES;
     DATAROW=2; 
RUN;

data open4;
set fake_open_data;
run;


proc format;
  value gendfmt
    1 = "Male"
    2 = "Female";
  value covfmt
    1 = "Var(Q)"
    2 = "Var(R)"
    3 = "Var(M)"
    4 = "Cov(Q,M)"
    5 = "Cov(R,M)"
    6 = "Cov(M1,M2)"
    7 = "Cov(Q,R)"
    8 = "Cov(Q1,Q2)"
    9 = "Cov(R1,R2)";
  run;

 /* get open data. */

data open;
  set open4;

  gender = gender + 1;   /* 1=male, 2=female. */

  log_bmi = log(bmi);

 /* log transform energy, protein and protein density. */
 /* m_ = biomarker, rec_ = recall, dhq_ = ffq.           */
 /* e = energy, p = protein, pd = protein density. */

  array x (*) m_e1 m_e2 rec_e1 rec_e2 dhq_e1 dhq_e2
              m_p1 m_p2 rec_p1 rec_p2 dhq_p1 dhq_p2
              m_pd1 m_pd2 rec_pd1 rec_pd2 dhq_pd1 dhq_pd2;

  do i = 1 to dim(x);
    if (x(i) ^= .) then x(i) = log(x(i));
    end;

  format gender gendfmt.;
  run;

proc sort data=open;  by gender id;  run;

 /* create data for proc mixed. */

data mixdata;
  set open;
  array yi (*) dhq_p1 dhq_p2 rec_p1 rec_p2 m_p1 m_p2;   /* ffq, recall and biomarker protein */

  /* create fixed covariates for mixed model. */

  array intercept (*) Intercept_Q1 Intercept_Q2 Intercept_R1 Intercept_R2 Intercept_M;
  array covar1 (*) age log_bmi;
  array covar2 (*) Age_Q Age_R Age_M
                   LogBMI_Q LogBMI_R LogBMI_M;

  obs = 0;

  do instr = 1 to 3;
    do repeat = 1 to 2;
      obs = obs + 1;
      y = yi(obs);

  /* create fixed covariates for mixed model. */

      do i = 1 to dim(intercept);  intercept(i) = 0;  end;
      if (instr <= 2) then intercept(obs) = 1;
      else intercept(5) = 1;

      do i = 1 to dim(covar2);  covar2(i) = 0;  end;
      do i = 1 to dim(covar1);
        covar2(3*(i-1) + instr) = covar1(i);
        end;

      output;
      end;
    end;
  run;

 /* specify covariance structure for proc mixed. */

 data ldata;
   parm = 1;
   row=1;  col1=1; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=0; col2=1; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   parm = 2;
   row=1;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=1; col4=0; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=0; col4=1; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   parm = 3;
   row=1;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=1; col6=0; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=1; output;
   parm = 4;
   row=1;  col1=0; col2=0; col3=0; col4=0; col5=1; col6=1; output;
   row=2;  col1=0; col2=0; col3=0; col4=0; col5=1; col6=1; output;
   row=3;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=5;  col1=1; col2=1; col3=0; col4=0; col5=0; col6=0; output;
   row=6;  col1=1; col2=1; col3=0; col4=0; col5=0; col6=0; output;
   parm = 5;
   row=1;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=0; col4=0; col5=1; col6=1; output;
   row=4;  col1=0; col2=0; col3=0; col4=0; col5=1; col6=1; output;
   row=5;  col1=0; col2=0; col3=1; col4=1; col5=0; col6=0; output;
   row=6;  col1=0; col2=0; col3=1; col4=1; col5=0; col6=0; output;
   parm = 6;
   row=1;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=1; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=1; col6=0; output;
   parm = 7;
   row=1;  col1=0; col2=0; col3=1; col4=1; col5=0; col6=0; output;
   row=2;  col1=0; col2=0; col3=1; col4=1; col5=0; col6=0; output;
   row=3;  col1=1; col2=1; col3=0; col4=0; col5=0; col6=0; output;
   row=4;  col1=1; col2=1; col3=0; col4=0; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   parm = 8;
   row=1;  col1=0; col2=1; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=1; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   parm = 9;
   row=1;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=2;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=3;  col1=0; col2=0; col3=0; col4=1; col5=0; col6=0; output;
   row=4;  col1=0; col2=0; col3=1; col4=0; col5=0; col6=0; output;
   row=5;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   row=6;  col1=0; col2=0; col3=0; col4=0; col5=0; col6=0; output;
   run;

 /* define covariate structure by gender. */

data ldata;
  set ldata;
  gender = 1;  output;
  gender = 2;  output;
  run;

proc sort data=ldata;  by gender parm row;  run;

 /* call proc mixed to fit model. */

proc mixed data=mixdata method=ml covtest cl asycov nobound;
  by gender;
  class obs;
  model y = Intercept_Q1 Intercept_Q2 Intercept_R1 Intercept_R2 Intercept_M
            Age_Q Age_R Age_M
            LogBMI_Q LogBMI_R LogBMI_M
            / noint solution;
  repeated obs / type=lin(9) ldata=ldata subject=id;
  parms 1 1 1 0 0 0 0 0 0;
  ods exclude modelinfo classlevels parmsearch iterhistory asycov lrt tests3;
  ods output covparms=covparms;
  ods output asycov=asycov;
  run;

 /* print covariance parameters. */

data covparms;
  set covparms;
    by gender;
  if (first.gender) then Parameter = 0;
  Parameter + 1;
  format Parameter covfmt.;
  run;

proc print data=covparms;
  by gender;
  id Parameter;
  var estimate stderr;
  title5 "Covariance Parameters";
  run;

title4;

 /* estimate attenuation factor and correlation with true intake. */

proc iml;
  gender = {"Male", "Female"};

  Var_T     = repeat(.,2,1);
  Corr_QT   = Var_T;
  Corr_RT   = Var_T;
  Corr_MT   = Var_T;
  Atten_Q   = Var_T;
  Atten_R   = Var_T;

  StdErr_VT = Var_T;
  StdErr_QT = Var_T;
  StdErr_RT = Var_T;
  StdErr_MT = Var_T;
  StdErr_AQ = Var_T;
  StdErr_AR = Var_T;

  do gend = 1 to 2;

    use covparms;
      read all var{estimate} where(gender=gend) into param;
      close covparms;

    use asycov;
      read all var("covp1":"covp6") where (gender=gend) into covar;
      close asycov;

    covar = covar[1:ncol(covar),];

    varQ   = param[1];
    varR   = param[2];
    varM   = param[3];
    covQM  = param[4];
    covRM  = param[5];
    varT   = param[6];

    Var_T[gend] = varT;
    StdErr_VT[gend] = sqrt(covar[6,6]);

    c1 = covQM / sqrt(varQ # varT);
	d1 = -.5 # (c1 / varQ) || {0 0} || (c1 / covQM) || {0} || -.5 # (c1 / varT);
    var1 = d1 * covar * d1`;

    Corr_QT[gend]   = c1;
	StdErr_QT[gend] = sqrt(var1);

    c2 = covRM / sqrt(varR # varT);
	d2 = {0} || -.5 # (c2 / varR) || {0 0} || (c2 / covRM) || -.5 # (c2 / varT);
    var2 = d2 * covar * d2`;

    Corr_RT[gend]   = c2;
	StdErr_RT[gend] = sqrt(var2);

    c3 = sqrt(varT / varM);
	d3 = {0 0} || -.5 # (c3 / varM) || {0 0} || .5 # (c3 / varT);
    var3 = d3 * covar * d3`;

    Corr_MT[gend]   = c3;
	StdErr_MT[gend] = sqrt(var3);

    a1 = covQM / varQ;
	d1 = -(a1 / varQ) || {0 0} || (a1 / covQM) || {0 0};
    var1 = d1 * covar * d1`;

    Atten_Q[gend]   = a1;
	StdErr_AQ[gend] = sqrt(var1);

    a2 = covRM / varR;
	d2 = {0} || -(a2 / varR) || {0 0} || (a2 / covRM) || {0};
    var2 = d2 * covar * d2`;

    Atten_R[gend]   = a2;
	StdErr_AR[gend] = sqrt(var2);
    end;

  print "Variance of True Intake",,
        gender Var_T StdErr_VT;
  print "Correlation of FFQ and True Intake",,
        gender Corr_QT stderr_QT;
  print "Correlation of 24HR and True Intake",,
        gender Corr_RT StdErr_RT;
  print "Correlation of UN and True Intake",,
        gender Corr_MT StdErr_MT;
  print "Attenuation Factor for FFQ",,
        gender Atten_Q StdErr_AQ;
  print "Attenuation Factor for 24HR",,
        gender Atten_R StdErr_AR;
  quit;