%macro SNP_NLMM(dataset = , subject = , response = , prog_stat = , log_cond_dens = , start_value = , q = , Kmax = 2, all_K = 1, info_crit = "BIC", quad_pt_K0 = 0, quad_pt_K1 = 0, grid_pt_K1 = 9, max_K1 = 1, quad_pt_K2 = 0, grid_pt_K2 = 9, max_K2 = 1, quad_pt_K3 = 0, grid_pt_K3 = 9, max_K3 = 1, quad_pt_K4 = 0, grid_pt_K4 = 9, max_K4 = 1); ods listing close; proc sort data=&dataset; by &subject; run; data &dataset; set &dataset; by &subject; last_obs=last.&subject; run; /**************************************************************/ /* Code for q=1 */ /**************************************************************/ %if &q=1 %then %do; /*************/ /* K=0 analysis using NLMIXED*/ /*************/ title "K=0 Analysis "; run; proc nlmixed data=&dataset %if &quad_pt_K0>0 %then %do; qpoints=&quad_pt_K0 %end; ; pi = 2*arsin(1); parms &start_value; /*** Define log conditional density of the response given the random effects ***/ /* b_i=mu+RZ_i */ b_1 = mu_1+(r_11*z_1); /* Additional programming statements */ &prog_stat; loglike_contribution=&log_cond_dens; model &response~general(loglike_contribution); random z_1 ~ normal([0],[1]) subject=&subject; predict r_11*z_1 out=z1_blup; estimate "expect_b1" mu_1; estimate "var_b1 " r_11**2; ods output ParameterEstimates=parm_k0 AdditionalEstimates=ae_k0 FitStatistics=fitstat_k0 Dimensions=dim_k0; run; /* Transpose parameter estimates into a dataset with one observation and columns for each variable */ proc transpose data=parm_k0 out=tparm_k0; var Estimate; id Parameter; run; proc transpose data=ae_k0 out=tae_k0; var Estimate; id Label; run; data tparm_k0; merge tparm_k0 tae_k0; test=1; run; /* Create dataset parm_total that will store all the parameter estimates for each K */ proc sort data=parm_k0 out=parm_total; by Parameter; run; data parm_total; set parm_total (keep=Parameter Estimate StandardError); rename Estimate=Est_K0; rename StandardError=StdErr_K0; run; /* Create dataset fitstat_total that will store all the information criteria and fit statistics for each K */ data fitstat_total; set fitstat_k0; rename Value=Value_K0; run; /* Create dataset dim_total that will store all the dimension information for each K */ data dim_total; set dim_k0; rename Value=Value_K0; run; /* Create dataset aw_total that will store all the additional estimates for each K */ proc sort data=ae_k0 out=ae_total; by Label; run; data ae_total; set ae_total (keep=Label Estimate StandardError); rename Estimate=Est_K0; rename StandardError=StdErr_K0; run; /* Evaluate the estimated density for the random effect - needed to plot the data */ /* Store the results in dataset graphics */ proc sql noprint; select expect_b1 into: graph_center from tparm_k0; quit; proc sql noprint; select var_b1 into: graph_spread from tparm_k0; quit; data graphics; do x_value = (&graph_center-4*sqrt(&graph_spread)) to (&graph_center+4*sqrt(&graph_spread)) by (sqrt(&graph_spread)*0.01); output; end; run; data graphics; set graphics; test=1; run; data graphics; merge graphics tparm_k0; by test; run; data graphics; set graphics; pi = 2*arsin(1); density_K0 = 1/sqrt(2*pi*var_b1)*exp(-1/(2*var_b1)*(x_value-expect_b1)**2); run; data graphics; set graphics (keep=x_value density_K0 test); run; /* For K>0 we will center and scale the numerical integration using the emprical Bayes estimates and corresponding standard error from assuming Gaussian random effects*/ /* Here we append the BLUPS and standard errors of b_1 assuming Gaussian random effects to the original dataset*/ data subjspec; set z1_blup (keep=&subject Pred StdErrPred last_obs); where last_obs=1; test=1; rename Pred=z1_blup; z1_var=StdErrPred**2; drop StdErrPred; run; data subjspec; merge subjspec tparm_k0; by test; G11 = sqrt(z1_var); b_hat_1 = z1_blup+expect_b1; run; data subjspec; set subjspec (keep=&subject G11 b_hat_1); run; data &dataset._2; merge &dataset subjspec; by &subject; run; /*************/ /* K=r analysis */ /* Note that we will loop over the values of K */ /*************/ %if &all_k=1 %then %do; %let start=1; %end; %if &all_k=0 %then %do; %let start=&Kmax; %end; %do r=&start %to &Kmax; /**************/ /* Here we give several items that are defined for each value of K including: 1) Grid of values of the vector xi 2) Number of quadrature points used in the numberical integration 3) Number of starting values that are used in the maximization 4) Bounds for the elements of xi 5) Construct c from the polar coordinate transformation for each of the grid points of xi 6) Construct a=B.inverse*c for each value of the grid points 7) Code for P_K() 8) Code to generate the expected value and variance of z_1 */ /**************/ %if &r=1 %then %do; /* Step 1 */ data xi_grid_q1; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K1-1)); output; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K1; /* Step 3 */ %let max_K=&max_K1; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1=sin(xi_1); c2=cos(xi_1);); /* Step 6 */ %let a_vector = %str(a0=c1; a1=c2;); /* Step 7 */ %let p_k = %str(p_k=(a0+a1*z_1);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*a0*a1; expect_z1sq = a0**2+3*a1**2; var_z1 =(expect_z1sq-expect_z1**2);); %end; %if &r=2 %then %do; /* Step 1 */ data xi_grid_q1; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); do xi_2 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); output; end; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K2; /* Step 3 */ %let max_K=&max_K2; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707 , -1.5707<= xi_2 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1 = sin(xi_1); c2 = cos(xi_1)*sin(xi_2); c3 = cos(xi_1)*cos(xi_2);); /* Step 6 */ %let a_vector = %str(a0=1.1944776*c1+-0.2705981*c3; a1=c2; a2=-0.2705981*c1+0.6532815*c3;); /* Step 7 */ %let p_k = %str(p_k=(a0+a1*z_1+a2*z_1**2);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*a0*a1+6*a1*a2; expect_z1sq = a0**2+3*a1**2+15*a2**2+6*a0*a2; var_z1=(expect_z1sq-expect_z1**2);); %end; %if &r=3 %then %do; /* Step 1 */ data xi_grid_q1; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_2 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_3 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); output; end; end; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K3; /* Step 3 */ %let max_K=&max_K3; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707 , -1.5707<= xi_2 <= 1.5707, -1.5707<= xi_3 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1 = sin(xi_1); c2 = cos(xi_1)*sin(xi_2); c3 = cos(xi_1)*cos(xi_2)*sin(xi_3); c4 = cos(xi_1)*cos(xi_2)*cos(xi_3);); /* Step 6 */ %let a_vector = %str(a0=1.1944776*c1+-0.2705981*c3; a1=1.5582767*c2+-0.2679064*c4; a2=-0.2705981*c1+0.6532815*c3; a3=-0.2679064*c2+0.3080468*c4;); /* Step 7 */ %let p_k = %str(p_k=(a0+a1*z_1+a2*z_1**2+a3*z_1**3);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*a0*a1+6*a1*a2+6*a0*a3+30*a2*a3; expect_z1sq = a0**2+3*a1**2+15*a2**2+105*a3**2+6*a0*a2+30*a1*a3; var_z1=(expect_z1sq-expect_z1**2);); %end; %if &r=4 %then %do; /* Step 1 */ data xi_grid_q1; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K4-1)); do xi_2 = -1.5 to 1.5 by (3/(&grid_pt_K4-1)); do xi_3 = -1.5 to 1.5 by (3/(&grid_pt_K4-1)); do xi_4 = -1.5 to 1.5 by (3/(&grid_pt_K4-1)); output; end; end; end; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K4; /* Step 3 */ %let max_K=&max_K4; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707 , -1.5707<= xi_2 <= 1.5707, -1.5707<= xi_3 <= 1.5707, -1.5707<= xi_4 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1 = sin(xi_1); c2 = cos(xi_1)*sin(xi_2); c3 = cos(xi_1)*cos(xi_2)*sin(xi_3); c4 = cos(xi_1)*cos(xi_2)*cos(xi_3)*sin(xi_4); c5 = cos(xi_1)*cos(xi_2)*cos(xi_3)*cos(xi_4);); /* Step 6 */ %let a_vector = %str(a0=1.28273545*c1+-0.4779613*c3+0.03380517*c5; a1=1.5582767*c2+-0.2679064*c4; a2=-0.47796128*c1+1.3210538*c3+-0.16238779*c5; a3=-0.2679064*c2+ 0.3080468*c4; a4=0.03380517*c1+ -0.1623878*c3+0.11897093*c5;); /* Step 7 */ %let p_k = %str(p_k=(a0+a1*z_1+a2*z_1**2+a3*z_1**3+a4*z_1**4);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*a0*a1+6*a1*a2+6*a0*a3+30*a2*a3+30*a1*a4+210*a3*a4; expect_z1sq = a0**2+3*a1**2+15*a2**2+105*a3**2+945*a4**2+6*a0*a2+30*a1*a3+30*a0*a4+210*a2*a4; var_z1=(expect_z1sq-expect_z1**2);); %end; /* Create data set of starting values where starting values for fixed-effect coefficients are taken from the estimates from assuming Gaussian random effects and mu_1 and r_11 are set so that the value of expect_b1 and var_b1 are the same as the estimates from assuming Gaussian random effects */ data xi_grid_q1; set xi_grid_q1; test=1; run; data start_value; merge tparm_k0 xi_grid_q1; by test; run; data start_value; set start_value; /* Construct c from the polar coordinate transformation */ &c_vector; /* Construct a=B.inverse*c for each value of the grid points */ &a_vector; /* Code to generate the expected value and variance of z_1 */ ℰ r_11 = sqrt(var_b1/var_z1); mu_1 = expect_b1-r_11*expect_z1; run; /* Drop values of a and c */ data start_value; set start_value; drop expect_z1 expect_z1sq var_z1 expect_b1 var_b1 test _NAME_ %do i=1 %to (&r+1); c&i %end; %do i=0 %to (&r); a&i %end; ; run; /*************/ /* K=r analysis using NLMIXED /*************/ /* First evaluate loglikelihood at each starting value */ /* No optimization is performed here */ /* Output the loglikelihood at each starting value in the dataset sv_k&r */ title "K=&r Analysis "; title2 "Log likelihood at starting values"; run; proc nlmixed data=&dataset._2 noad tech=none %if &quad_pt >0 %then %do; qpoints=&quad_pt %end; ; pi = 2*arsin(1); parms /data=start_value; bounds &bounds_K; /* Z_i=R.inv(b_hat-mu+G^{1/2}U) */ J = b_hat_1-mu_1+G11*u_1; z_1 = J/r_11; /* Construct c from the polar coordinate transformation */ &c_vector; /* Construct a=B.inverse*c for each value of the grid points */ &a_vector; /* Define p_k in the RE density */ &p_k; /* Define the log random effects density */ log_re_dens = log((p_k)**2)-0.5*(z_1**2)-1/2*log(2*pi); /* Determinant of square root of G matrix */ log_det_G = log((G11)); /* Determinant of inverse of R */ log_det_R_inv = log(1/(r_11)); /* Define random effects density with respect to u_i */ log_re_dens = log_re_dens+log_det_G+log_det_R_inv; /*** Define log conditional density of the response given random effects ***/ /* b_i=mu+RZ_i */ b_1 = mu_1+(r_11*z_1); /* Additional programming statements */ &prog_stat; /*** Give the log likelihood reformulation ***/ /* log(random effects density) - log(standard q-dimensional normal) */ log_re_reform = log_re_dens-(-0.5*(u_1**2)-1/2*log(2*pi)); /* Give log likelihood contribution */ if last_obs=0 then loglike_contribution=&log_cond_dens; if last_obs=1 then loglike_contribution=&log_cond_dens+log_re_reform; model &response~general(loglike_contribution); random u_1 ~ normal([0],[1]) subject=&subject; ods output Parameters=sv_k&r; run; /* Sort sv_k&r by log likelihood */ proc sort data=sv_k&r; by NegLogLike; run; data sv_k&r; set sv_k&r; rank = _N_; run; /*********************/ /* Maximize the likelihood for each of the set of staring values with the max_k^{th} smallest log likelihood */ /*********************/ %do i=1 %to &max_k; run; data sv_k&r._&i; set sv_k&r; where rank=&i; run; data sv_k&r._&i; set sv_k&r._&i; drop NegLogLike rank; run; title "K=&r Analysis "; title2 "Maximization using starting value set &i"; run; proc nlmixed data=&dataset._2 noad %if &quad_pt >0 %then %do; qpoints=&quad_pt %end; ; pi = 2*arsin(1); parms /data=sv_k&r._&i; bounds &bounds_K; /* Z_i=R.inv(b_hat-mu+G^{1/2}U) */ J = b_hat_1-mu_1+G11*u_1; z_1 = J/r_11; /* Construct c from the polar coordinate transformation */ &c_vector; /* Construct a=B.inverse*c for each value of the grid points */ &a_vector; /* Define p_k in the RE density */ &p_k; /* Define the log random effects density */ log_re_dens = log((p_k)**2)-0.5*(z_1**2)-1/2*log(2*pi); /* Determinant of square root of G matrix */ log_det_G = log((G11)); /* Determinant of inverse of R */ log_det_R_inv = log(1/(r_11)); /* Define random effects density with respect to u_i */ log_re_dens = log_re_dens+log_det_G+log_det_R_inv; /*** Define log conditional density of the response given random effects ***/ /* b_i=mu+RZ_i */ b_1 = mu_1+(r_11*z_1); /* Additional programming statements */ &prog_stat; /*** Give the log likelihood reformulation ***/ /* i.e. log(random effects density) - log(standard q-dimensional normal) */ log_re_reform = log_re_dens-(-0.5*(u_1**2)-1/2*log(2*pi)); /* Give log likelihood contribution */ if last_obs=0 then loglike_contribution=&log_cond_dens; if last_obs=1 then loglike_contribution=&log_cond_dens+log_re_reform; model &response~general(loglike_contribution); random u_1 ~ normal([0],[1]) subject=&subject; /*** Estimate statements needed to obtain the mean and variance of b_i ***/ /* Compute the mean and variance of Z_i for each K */ ℰ /* E(b_i)=mu+R*E(Z_i) */ estimate "expect_b1" mu_1+r_11*expect_z1; /* var(b_i)=R*var(Z_i)*R^T */ estimate "var_b1" r_11**2*var_z1; ods output ParameterEstimates=parm_k&r._&i AdditionalEstimates=ae_k&r._&i FitStatistics=fitstat_k&r._&i Dimensions=dim_k&r._&i; run; %end; /* End loop over the starting values */ /* Choose the final parameter estimates which give the smallest log likelihood */ data loglike_k&r; set %do i=1 %to &max_k; fitstat_k&r._&i %end; ; run; data loglike_k&r; set loglike_k&r; where Descr="-2 Log Likelihood"; q_iteration = _N_; run; proc sort data=loglike_k&r; by Value; run; data loglike_k&r; set loglike_k&r; rank = _N_; data loglike_min_k&r; set loglike_k&r; where rank=1; run; proc sql noprint; select q_iteration into: q_it from loglike_min_k&r; quit; %do i=&q_it %to &q_it; data parm_k&r; set parm_k&r._&i; run; data ae_k&r; set ae_k&r._&i; run; data fitstat_k&r; set fitstat_k&r._&i; run; data dim_k&r; set dim_k&r._&i; run; %end; /* Transpose parameter estimates into a dataset with one observation and columns for each variable */ proc transpose data=parm_k&r out=tparm_k&r; var Estimate; id Parameter; run; proc transpose data=ae_k&r out=tae_k&r; var Estimate; id Label; run; data tparm_k&r; merge tparm_k&r tae_k&r; test=1; run; /* Evaluate the estimated density for the random effect - needed to plot the data */ /* Store the results in dataset graphics_k&r */ data graphics_k&r; merge graphics tparm_k&r; by test; run; data graphics_k&r; set graphics_k&r; pi = 2*arsin(1); z_1 = (x_value-mu_1)/r_11; /* Step 1 */ &c_vector; /* Step 2 */ &a_vector; /* Step 3 */ &p_k; density_K&r = (p_k**2)*1/sqrt(2*pi)*exp(-1/2*(z_1)**2)/r_11; run; data graphics_K&r; set graphics_K&r (keep=x_value density_K&r test); run; /* Add to dataset parm_total with parameter estimates from each K=r */ proc sort data=parm_k&r; by Parameter; run; data sparm_k&r; set parm_k&r (keep=Parameter Estimate StandardError); rename Estimate=Est_K&r; rename StandardError=StdErr_K&r; run; data parm_total; merge parm_total sparm_k&r; by Parameter; run; /* Add to dataset ae_total with additional estimates from each K=r */ proc sort data=ae_k&r; by Label; run; data sae_k&r; set ae_k&r (keep=Label Estimate StandardError); rename Estimate=Est_K&r; rename StandardError=StdErr_K&r; run; data ae_total; merge ae_total sae_k&r; by Label; run; /* Add to dataset fitstat_total with information criteria and fit statistics from each K=r */ data fitstat_k&r; set fitstat_k&r; rename Value=Value_K&r; run; data fitstat_total; merge fitstat_total fitstat_k&r; by Descr; run; data dim_k&r; set dim_k&r; rename Value=Value_K&r; run; /* Add to dataset dim_total with dimensions from each K=r */ data dim_total; merge dim_total dim_k&r; run; %end; /* End looping over values of K */ /****************/ /* Print Results of Interest */ /****************/ ods listing; data ae_total; set ae_total; rename Label=Parameter; run; data parm_total; set ae_total parm_total ; run; data fitstat_total; set fitstat_total; rename Descr=Criteria; run; /* Determine which K is selected by the preferred information criteria */ %if &info_crit="BIC" %then %do; data fit_crit; set fitstat_total; where Criteria="BIC (smaller is better)"; run; %end; %if &info_crit="AIC" %then %do; data fit_crit; set fitstat_total; where Criteria="AIC (smaller is better)"; run; %end; proc transpose data=fit_crit out=tfit_crit; by criteria; run; data tfit_crit; set tfit_crit; K_number=(_N_ -1); run; proc sort data=tfit_crit; by COL1; run; data tfit_crit; set tfit_crit; rank = _N_; run; data tfit_crit; set tfit_crit; where rank=1; run; proc sql noprint; select K_number into: K_select from tfit_crit; quit; %if &all_K=0 & &K_select=1 %then %do; %let K_select=&Kmax; %end; /* For the K selected combine the parameter estimates and additional estimates in one dataset */ %do r=&K_select %to &K_select; data ae_k&r; set ae_k&r; rename Label=Parameter; run; data parm_final_K; set ae_k&r parm_k&r; run; %end; title "Model Information"; title2 "By Degree of Flexibility"; proc print data=dim_total noobs; run; title "Information Criteria"; title2 "By Degree of Flexibility"; proc print data=fitstat_total noobs; run; title "Parameter Estimates"; title2 "By Degree of Flexibility"; proc print data=parm_total noobs; run; title "Parameter Estimates"; title2 "K = &K_select"; title3 "K selected by"; title4 &info_crit; proc print data=parm_final_K noobs; run; /*******************/ /* Graph estimated densities */ /*******************/ symbol1 interpol=join value=none line=1 CI=blue; axis1 label=("Random Coefficient"); axis2 label=("Density"); title "K=0 Estimated Random Effects Density"; proc gplot data=graphics; plot density_K0*x_value / haxis=axis1 vaxis=axis2 ; run; %do r=(1*&all_K+&Kmax*(1-&all_K)) %to &Kmax; title "K=&r Estimated Random Effects Density"; proc gplot data=graphics_K&r; plot density_K&r *x_value / haxis=axis1 vaxis=axis2 ; run; %end; quit; %end; /*End q=1 Code*/ /**************************************************************/ /* Code for q=2 */ /**************************************************************/ %if &q=2 %then %do; /*************/ /* K=0 analysis using NLMIXED*/ /*************/ title "K=0 Analysis "; run; proc nlmixed data=&dataset %if &quad_pt_K0>0 %then %do; qpoints=&quad_pt_K0 %end; ; pi = 2*arsin(1); parms &start_value; /*** Define log conditional density of response given random effects ***/ b_1 = mu_1+(r_11*z_1); b_2 = mu_2+(r_21*z_1+r_22*z_2); /* Additional programming statements */ &prog_stat; loglike_contribution=&log_cond_dens; model &response~general(loglike_contribution); random z_1 z_2 ~ normal([0,0],[1,0,1]) subject=&subject; predict r_11*z_1 out=z1_blup; predict r_21*z_1+r_22*z_2 out=z2_blup; predict r_11*z_1+r_21*z_1+r_22*z_2 out=z1z2_blup; estimate "expect_b1" mu_1; estimate "expect_b2" mu_2; estimate "var_b1 " r_11**2; estimate "cov_b1b2" r_11*r_21; estimate "var_b2" r_21**2+r_22**2; ods output ParameterEstimates=parm_k0 AdditionalEstimates=ae_k0 FitStatistics=fitstat_k0 Dimensions=dim_k0; run; /* Transpose parameter estimates into a dataset with one observation and columns for each variable */ proc transpose data=parm_k0 out=tparm_k0; var Estimate; id Parameter; run; proc transpose data=ae_k0 out=tae_k0; var Estimate; id Label; run; data tparm_k0; merge tparm_k0 tae_k0; test=1; run; /* Create dataset parm_total that will store all the parameter estimates for each K */ proc sort data=parm_k0 out=parm_total; by Parameter; run; data parm_total; set parm_total (keep=Parameter Estimate StandardError); rename Estimate=Est_K0; rename StandardError=StdErr_K0; run; /* Create dataset fitstat_total that will store all the information criteria and fit statistics for each K */ data fitstat_total; set fitstat_k0; rename Value=Value_K0; run; /* Create dataset dim_total that will store all the dimension information for each K */ data dim_total; set dim_k0; rename Value=Value_K0; run; /* Create dataset aw_total that will store all the additional estimates for each K */ proc sort data=ae_k0 out=ae_total; by Label; run; data ae_total; set ae_total (keep=Label Estimate StandardError); rename Estimate=Est_K0; rename StandardError=StdErr_K0; run; /* Evaluate the estimated density for the random effect - needed to plot the data */ /* Store the results in dataset graphics */ proc sql noprint; select expect_b1 into: graph_center1 from tparm_k0; quit; proc sql noprint; select var_b1 into: graph_spread1 from tparm_k0; quit; proc sql noprint; select expect_b2 into: graph_center2 from tparm_k0; quit; proc sql noprint; select var_b2 into: graph_spread2 from tparm_k0; quit; data graphics; do x_value1 = (&graph_center1-4*sqrt(&graph_spread1)) to (&graph_center1+4*sqrt(&graph_spread1)) by (sqrt(&graph_spread1)*0.05); do x_value2 = (&graph_center2-4*sqrt(&graph_spread2)) to (&graph_center2+4*sqrt(&graph_spread2)) by (sqrt(&graph_spread2)*0.05); output; end; end; run; data graphics; set graphics; test=1; run; data graphics; merge graphics tparm_k0; by test; run; data graphics; set graphics; pi = 2*arsin(1); rho = cov_b1b2/sqrt(var_b1*var_b2); z = (x_value1-expect_b1)**2/var_b1-2*rho*(x_value1-expect_b1)*(x_value2-expect_b2)/sqrt(var_b1*var_b2)+ (x_value2-expect_b2)**2/var_b2; density_K0 = 1/(2*pi*sqrt(var_b1*var_b2*(1-rho**2)))*exp(-1/(2*(1-rho**2))*z); run; data graphics; set graphics (keep=x_value1 x_value2 density_K0 test); run; /* For K>0 we will center and scale the numerical integration using the emprical Bayes estimates and corresponding standard error from assuming Gaussian random effects*/ /* Here we append the BLUPS and standard errors of b_1 and b_2 assuming Gaussian random effects to the original dataset*/ data z1_blup; set z1_blup (keep=&subject Pred StdErrPred last_obs); where last_obs=1; rename Pred= z1_blup; z1_var =StdErrPred**2; drop StdErrPred; run; data z2_blup; set z2_blup (keep=&subject Pred StdErrPred last_obs); where last_obs=1; rename Pred= z2_blup; z2_var =StdErrPred**2; drop StdErrPred; run; data z1z2_blup; set z1z2_blup (keep=&subject Pred StdErrPred last_obs); where last_obs=1; rename Pred= z1z2_blup; z1z2_var =StdErrPred**2; drop StdErrPred; run; data subjspec; merge z1_blup z2_blup z1z2_blup; by &subject; test = 1; z1z2_cov = (z1z2_var-z1_var-z2_var)/2; drop last_obs; run; data subjspec; merge subjspec tparm_k0; by test; run; data subjspec; set subjspec; lambda1 = ((z1_var+z2_var)+sqrt((z2_var-z1_var)**2+4*z1z2_cov**2))/2; lambda2 = ((z1_var+z2_var)-sqrt((z2_var-z1_var)**2+4*z1z2_cov**2))/2; Q11 = sqrt(1/(1+(z1_var-lambda1)**2/z1z2_cov**2)); Q21 = -(z1_var-lambda1)*Q11/z1z2_cov; Q12 = sqrt(1/(1+(z1_var-lambda2)**2/z1z2_cov**2)); Q22 = -(z1_var-lambda2)*Q12/z1z2_cov; /* Note: Gij is the (i,j)th element of the square root matrix */ G11 = sqrt(lambda1)*Q11**2+sqrt(lambda2)*Q12*Q21; G12 = sqrt(lambda1)*Q11*Q12+sqrt(lambda2)*Q12*Q22; G21 = sqrt(lambda1)*Q11*Q21+sqrt(lambda2)*Q21*Q22; G22 = sqrt(lambda1)*Q12*Q21+sqrt(lambda2)*Q22**2; b_hat_1 = z1_blup+expect_b1; b_hat_2 = z2_blup+expect_b2; run; data subjspec; set subjspec (keep=&subject G11 G12 G21 G22 b_hat_1 b_hat_2); run; data &dataset._2; merge &dataset subjspec; by &subject; run; /*************/ /* K=r analysis */ /* Note that we will loop over the values of K */ /*************/ %if &all_k=1 %then %do; %let start=1; %end; %if &all_k=0 %then %do; %let start=&Kmax; %end; %do r=&start %to &Kmax; /**************/ /* Here we give several items that are defined for each value of K including: 1) Grid of values of the vector xi 2) Number of quadrature points used in the numberical integration 3) Number of starting values that are used in the maximization 4) Bounds for the elements of xi 5) Construct c from the polar coordinate transformation for each of the grid points of xi 6) Construct a=B.inverse*c for each value of the grid points 7) Code for P_K() 8) Code to generate the expected value and variance of z_1 */ /**************/ %if &r=1 %then %do; /* Step 1 */ data xi_grid_q2; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K1-1)); do xi_2 = -1.5 to 1.5 by (3/(&grid_pt_K1-1)); output; end; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K1; /* Step 3 */ %let max_K=&max_K1; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707, -1.5707<= xi_2 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1=sin(xi_1); c2=cos(xi_1)*sin(xi_2); c3=cos(xi_1)*cos(xi_2);); /* Step 6 */ %let a_vector = %str(a00=c1; a10=c2; a01=c3;); /* Step 7 */ %let p_k = %str(p_k=(a00+a10*z_1+a01*z_2);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*a00*a10; expect_z2 = 2*a00*a01; expect_z1sq = a00**2+3*a10**2+a01**2; expect_z1z2 = 2*a10*a01; expect_z2sq = a00**2+a10**2+3*a01**2; var_z1 =(expect_z1sq-expect_z1**2); cov_z1z2 =(expect_z1z2-expect_z1*expect_z2); var_z2 =(expect_z2sq-expect_z2**2);); %end; %if &r=2 %then %do; /* Step 1 */ data xi_grid_q2; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); do xi_2 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); do xi_3 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); do xi_4 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); do xi_5 = -1.5 to 1.5 by (3/(&grid_pt_K2-1)); output; end; end; end; end; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K2; /* Step 3 */ %let max_K=&max_K2; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707 , -1.5707<= xi_2 <= 1.5707, -1.5707<= xi_3 <= 1.5707, -1.5707<= xi_4 <= 1.5707, -1.5707<= xi_5 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1=sin(xi_1); c2=cos(xi_1)*sin(xi_2); c3=cos(xi_1)*cos(xi_2)*sin(xi_3); c4=cos(xi_1)*cos(xi_2)*cos(xi_3)*sin(xi_4); c5=cos(xi_1)*cos(xi_2)*cos(xi_3)*cos(xi_4)*sin(xi_5); c6=cos(xi_1)*cos(xi_2)*cos(xi_3)*cos(xi_4)*cos(xi_5);); /* Step 6 */ %let a_vector = %str(a00=1.3683057*c1-0.25272473*c4-0.25272473*c6; a10=c2; a01=c3; a20=-0.25272473*c1+0.65861913*c4-0.04848765*c6; a11=c5; a02=-0.25272473*c1+0.65861913*c6-0.04848765*c4;); /* Step 7 */ %let p_k = %str(p_k=(a00+a10*z_1+a01*z_2+a20*z_1**2+a11*z_1*z_2+a02*z_2**2);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*a00*a10+2*a11*a01+2*a10*a02+6*a20*a10; expect_z2 = 2*a00*a01+2*a10*a11+2*a01*a20+6*a01*a02; expect_z1sq = a00**2+3*a10**2+a01**2+15*a20**2+3*a11**2+ 3*a02**2+6*a00*a20+6*a20*a02+2*a00*a02; expect_z1z2 = 2*a10*a01+2*a00*a11+6*a20*a11+6*a11*a02; expect_z2sq = a00**2+a10**2+3*a01**2+3*a20**2+3*a11**2+ 15*a02**2+2*a00*a20+6*a00*a02+6*a20*a02; var_z1 =(expect_z1sq-expect_z1**2); cov_z1z2 =(expect_z1z2-expect_z1*expect_z2); var_z2 =(expect_z2sq-expect_z2**2);); %end; %if &r=3 %then %do; /* Step 1 */ data xi_grid_q2; do xi_1 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_2 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_3 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_4 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_5 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_6 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_7 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_8 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); do xi_9 = -1.5 to 1.5 by (3/(&grid_pt_K3-1)); output; end; end; end; end; end; end; end; end; end; run; /* Step 2 */ %let quad_pt=&quad_pt_K3; /* Step 3 */ %let max_K=&max_K3; /* Step 4 */ %let bounds_K = %str(-1.5707<= xi_1 <= 1.5707 , -1.5707<= xi_2 <= 1.5707, -1.5707<= xi_3 <= 1.5707, -1.5707<= xi_4 <= 1.5707, -1.5707<= xi_5 <= 1.5707 -1.5707<= xi_6 <= 1.5707, -1.5707<= xi_7 <= 1.5707, -1.5707<= xi_8 <= 1.5707, -1.5707<= xi_9 <= 1.5707;); /* Step 5 */ %let c_vector = %str(c1=sin(phi_1); c2=cos(phi_1)*sin(phi_2); c3=cos(phi_1)*cos(phi_2)*sin(phi_3); c4=cos(phi_1)*cos(phi_2)*cos(phi_3)*sin(phi_4); c5=cos(phi_1)*cos(phi_2)*cos(phi_3)*cos(phi_4)*sin(phi_5); c6=cos(phi_1)*cos(phi_2)*cos(phi_3)*cos(phi_4)*cos(phi_5)*sin(phi_6); c7=cos(phi_1)*cos(phi_2)*cos(phi_3)*cos(phi_4)*cos(phi_5)*cos(phi_6)*sin(phi_7); c8=cos(phi_1)*cos(phi_2)*cos(phi_3)*cos(phi_4)*cos(phi_5)*cos(phi_6)*cos(phi_7)*sin(phi_8); c9=cos(phi_1)*cos(phi_2)*cos(phi_3)*cos(phi_4)*cos(phi_5)*cos(phi_6)*cos(phi_7)*cos(phi_8)*sin(phi_9); c10=cos(phi_1)*cos(phi_2)*cos(phi_3)*cos(phi_4)*cos(phi_5)*cos(phi_6)*cos(phi_7)*cos(phi_8)*cos(phi_9);); /* Step 6 */ %let a_vector = %str(a00=1.3683057*c1-0.25272473*c4-0.25272473*c6; a10=1.6994858*c2-0.25430903*c7-0.21696726*c9; a01=1.6994858*c3-0.25430903*c10-0.21696726*c8; a20=-0.25272473*c1+0.65861913*c4-0.04848765*c6; a11=c5; a02=-0.25272473*c1+0.65861913*c6-0.04848765*c4; a30 = -0.2543090*c2+0.31441401*c7+-0.05601265*c9; a21 = -0.2169673*c3+0.67066220*c8-0.05601265*c10; a12 = -0.2169673*c2+0.67066220*c9-0.05601265*c7; a03 = -0.2543090*c3+0.31441401*c10+-0.05601265*c8;); /* Step 7 */ %let p_k = %str(p_k=(a00+a10*z_1+a01*z_2+a20*z_1**2+a11*z_1*z_2+a02*z_2**2+a30*z_1**3+a21*z_1**2*z_2+a12*z_1*z_2**2+a03*z_2**3);); /* Step 8 */ %let expectation = %str(expect_z1 = 2*(a00*a10+3*a00*a30+a00*a12+3*a10*a20+a10*a02+a01*a11+15*a20*a30+3*a20*a12+3*a11*a21+3*a11*a03+3*a02*a30+3*a02*a12); expect_z2 = 2*(a00*a01+a00*a21+3*a00*a03+a10*a11+a01*a20+3*a01*a02+3*a20*a21+3*a20*a03+3*a11*a30+3*a11*a12+3*a02*a21+15*a02*a03); expect_z1sq = a00**2+3*a10**2+a01**2+15*a20**2+3*a11**2+3*a02**2+105*a30**2+15*a21**2+9*a12**2+15*a03**2+ 2*(3*a00*a20+a00*a02+15*a10*a30+3*a10*a12+3*a01*a21+3*a01*a03+3*a20*a02+15*a30*a12+9*a21*a03); expect_z1z2 = 2*(a10*a01+a00*a11+3*a10*a21+3*a10*a03+3*a01*a30+3*a01*a12+3*a20*a11+3*a11*a02+15*a30*a21+9*a30*a03+9*a21*a12+15*a12*a03); expect_z2sq = a00**2+a10**2+3*a01**2+3*a20**2+3*a11**2+15*a02**2+15*a30**2+9*a21**2+15*a12**2+105*a03**2+ 2*(a00*a20+3*a00*a02+3*a10*a30+3*a10*a12+3*a01*a21+15*a01*a03+3*a20*a02+9*a30*a12+15*a21*a03); var_z1 =(expect_z1sq-expect_z1**2); cov_z1z2 =(expect_z1z2-expect_z1*expect_z2); var_z2 =(expect_z2sq-expect_z2**2);); %end; /* Create data set of starting values where starting values for fixed-effect coefficients are taken from the estimates from assuming Gaussian random effects and mu_1 and r_11 are set so that the value of expect_b1 and var_b1 are the same as the estimates from assuming Gaussian random effects */ data xi_grid_q2; set xi_grid_q2; test=1; run; data start_value; merge tparm_k0 xi_grid_q2; by test; run; data start_value; set start_value; /* Construct c from the polar coordinate transformation */ &c_vector; /* Construct a=B.inverse*c for each value of the grid points */ &a_vector; /* Code to generate the expected value and variance of z_1 */ ℰ r_11 = sqrt(var_b1/var_z1); B_q = 2*cov_z1z2*cov_b1b2/(var_z1*r_11)-2*cov_b1b2* cov_z1z2/( var_z1*r_11); A_q = (cov_z1z2/ var_z1)**2-2* cov_z1z2**2/ var_z1+var_z2; C_q = var_z1*(cov_b1b2/( var_z1*r_11))**2-var_b2; r_22 = (-B_q+sqrt(B_q**2-4*A_q*C_q))/(2*A_q); r_21 = (cov_b1b2-cov_z1z2*r_11*r_22)/(var_z1*r_11); /* mu = E(b_i)-RE(Z_i) */ mu_1 = expect_b1-r_11*expect_z1; mu_2 = expect_b2-r_21*expect_z1-r_22*expect_z2; run; /* Drop values of a and c */ data start_value; set start_value; drop expect_z1 expect_z2 expect_z1sq expect_z1z2 expect_z2sq var_z1 cov_z1z2 var_z2 expect_b1 expect_b2 var_b1 var_b2 cov_b1b2 B_q A_q C_q test _NAME_; run; %if &r=1 %then %do; data start_value; set start_value; drop c1 c2 c3 a00 a10 a01; run; %end; %if &r=2 %then %do; data start_value; set start_value; drop c1 c2 c3 c4 c5 c6 a00 a10 a01 a20 a11 a02; run; %end; %if &r=3 %then %do; data start_value; set start_value; drop c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 a00 a10 a01 a20 a11 a02 a30 a21 a12 a03; run; %end; /*************/ /* K=r analysis using NLMIXED /*************/ /* First evaluate loglikelihood at each starting value */ /* No optimization is performed here */ /* Output the loglikelihood at each starting value in the dataset sv_k&r */ title "K=&r Analysis "; title2 "Log likelihood at starting values"; run; proc nlmixed data=&dataset._2 noad tech=none %if &quad_pt >0 %then %do; qpoints=&quad_pt %end; ; pi = 2*arsin(1); parms /data=start_value; bounds &bounds_K; /* Z_i=R.inv(b_hat-mu+G^{1/2}U) */ J = b_hat_1-mu_1+G11*u_1+G12*u_2; K = b_hat_2-mu_2+G21*u_1+G22*u_2; z_1 = J/r_11; z_2 =-r_21/(r_11*r_22)*J+K/r_22; /* Construct c from the polar coordinate transformation */ &c_vector; /* Construct a=B.inverse*c for each value of the grid points */ &a_vector; /* Define p_k in the RE density */ &p_k; /* log SNP density = log(P_K^2)+log(standard q-dimensional normal) */ log_re_dens = log((p_k)**2)-0.5*(z_1**2+z_2**2)-2/2*log(2*pi); /* Determinant of square root of G matrix */ log_det_G = log((G11*G22-G12*G21)); /* Determinant of inverse of R */ log_det_R_inv = log(1/(r_11*r_22)); /* Define random effects density with respect to u_i */ log_re_dens = log_re_dens+log_det_G+log_det_R_inv; /*** Define log conditional density of the response given random effects ***/ /* b_i=mu+RZ_i */ b_1 = mu_1+(r_11*z_1); b_2 = mu_2+(r_21*z_1+r_22*z_2); /* Additional programming statements */ &prog_stat; /* log(random effects density) - log(standard q-dimensional normal) */ log_re_reform = log_re_dens-(-0.5*(u_1**2+u_2**2)-2/2*log(2*pi)); /* Give log likelihood contribution */ if last_obs=0 then loglike_contribution=&log_cond_dens; if last_obs=1 then loglike_contribution=&log_cond_dens+log_re_reform; model &response~general(loglike_contribution); random u_1 u_2 ~ normal([0,0],[1,0,1]) subject=&subject; ods output Parameters=sv_k&r; run; /* Sort sv_k&r by log likelihood */ proc sort data=sv_k&r; by NegLogLike; run; data sv_k&r; set sv_k&r; rank = _N_; run; /*********************/ /* Maximize the likelihood for each of the set of staring values with the max_k^{th} smallest log likelihood */ /*********************/ %do i=1 %to &max_k; run; data sv_k&r._&i; set sv_k&r; where rank=&i; run; data sv_k&r._&i; set sv_k&r._&i; drop NegLogLike rank; run; title "K=&r Analysis "; title2 "Maximization using starting value set &i"; run; proc nlmixed data=&dataset._2 noad %if &quad_pt >0 %then %do; qpoints=&quad_pt %end; ; pi = 2*arsin(1); parms /data=sv_k&r._&i; bounds &bounds_K; /* Z_i=R.inv(b_hat-mu+G^{1/2}U) */ J = b_hat_1-mu_1+G11*u_1+G12*u_2; K = b_hat_2-mu_2+G21*u_1+G22*u_2; z_1 = J/r_11; z_2 =-r_21/(r_11*r_22)*J+K/r_22; /* Construct c from the polar coordinate transformation */ &c_vector; /* Construct a=B.inverse*c for each value of the grid points */ &a_vector; /* Define p_k in the RE density */ &p_k; /* log SNP density = log(P_K^2)+log(standard q-dimensional normal) */ log_re_dens = log((p_k)**2)-0.5*(z_1**2+z_2**2)-2/2*log(2*pi); /* Determinant of square root of G matrix */ log_det_G = log((G11*G22-G12*G21)); /* Determinant of inverse of R */ log_det_R_inv = log(1/(r_11*r_22)); /* Define random effects density with respect to u_i */ log_re_dens = log_re_dens+log_det_G+log_det_R_inv; /*** Define log conditional density of the response given random effects ***/ /* b_i=mu+RZ_i */ b_1 = mu_1+(r_11*z_1); b_2 = mu_2+(r_21*z_1+r_22*z_2); /* Additional programming statements */ &prog_stat; /* log(random effects density) - log(standard q-dimensional normal) */ log_re_reform = log_re_dens-(-0.5*(u_1**2+u_2**2)-2/2*log(2*pi)); /* Give log likelihood contribution */ if last_obs=0 then loglike_contribution=&log_cond_dens; if last_obs=1 then loglike_contribution=&log_cond_dens+log_re_reform; model &response~general(loglike_contribution); random u_1 u_2 ~ normal([0,0],[1,0,1]) subject=&subject; /*** Estimate statements needed to obtain the mean and variance of b_i ***/ /* Compute the mean and variance of Z_i for each K */ ℰ /* E(b_i)=mu+R*E(Z_i) */ estimate "expect_b1" mu_1+r_11*expect_z1; estimate "expect_b2" mu_2+r_21*expect_z1+r_22*expect_z2; /* var(b_i)=R*var(Z_i)*R^T */ estimate "var_b1" r_11**2*var_z1; estimate "cov_b1b2" var_z1*r_11*r_21+cov_z1z2*r_11*r_22; estimate "var_b2" var_z1*r_21**2+2*cov_z1z2*r_21*r_22+var_z2*r_22**2; ods output ParameterEstimates=parm_k&r._&i AdditionalEstimates=ae_k&r._&i FitStatistics=fitstat_k&r._&i Dimensions=dim_k&r._&i; run; %end; /* End loop over the starting values */ /* Choose the final parameter estimates which give the smallest log likelihood */ data loglike_k&r; set %do i=1 %to &max_k; fitstat_k&r._&i %end; ; run; data loglike_k&r; set loglike_k&r; where Descr="-2 Log Likelihood"; q_iteration = _N_; run; proc sort data=loglike_k&r; by Value; run; data loglike_k&r; set loglike_k&r; rank = _N_; data loglike_min_k&r; set loglike_k&r; where rank=1; run; proc sql noprint; select q_iteration into: q_it from loglike_min_k&r; quit; %do i=&q_it %to &q_it; data parm_k&r; set parm_k&r._&i; run; data ae_k&r; set ae_k&r._&i; run; data fitstat_k&r; set fitstat_k&r._&i; run; data dim_k&r; set dim_k&r._&i; run; %end; /* Transpose parameter estimates into a dataset with one observation and columns for each variable */ proc transpose data=parm_k&r out=tparm_k&r; var Estimate; id Parameter; run; proc transpose data=ae_k&r out=tae_k&r; var Estimate; id Label; run; data tparm_k&r; merge tparm_k&r tae_k&r; test=1; run; /* Evaluate the estimated density for the random effect - needed to plot the data */ /* Store the results in dataset graphics_k&r */ data graphics_k&r; merge graphics tparm_k&r; by test; run; data graphics_k&r; set graphics_k&r; pi = 2*arsin(1); /* Z_i=R.inv(x_value-mu) */ J = x_value1-mu_1; K = x_value2-mu_2; z_1 = J/r_11; z_2 =-r_21/(r_11*r_22)*J+K/r_22; /* Step 1 */ &c_vector; /* Step 2 */ &a_vector; /* Step 3 */ &p_k; density_K&r = (p_k**2)*1/(2*pi)*exp(-1/2*(z_1**2+z_2**2))/abs(r_11*r_22); run; data graphics_K&r; set graphics_K&r (keep=x_value1 x_value2 density_K&r test); run; /* Add to dataset parm_total with parameter estimates from each K=r */ proc sort data=parm_k&r; by Parameter; run; data sparm_k&r; set parm_k&r (keep=Parameter Estimate StandardError); rename Estimate=Est_K&r; rename StandardError=StdErr_K&r; run; data parm_total; merge parm_total sparm_k&r; by Parameter; run; /* Add to dataset ae_total with additional estimates from each K=r */ proc sort data=ae_k&r; by Label; run; data sae_k&r; set ae_k&r (keep=Label Estimate StandardError); rename Estimate=Est_K&r; rename StandardError=StdErr_K&r; run; data ae_total; merge ae_total sae_k&r; by Label; run; /* Add to dataset fitstat_total with information criteria and fit statistics from each K=r */ data fitstat_k&r; set fitstat_k&r; rename Value=Value_K&r; run; data fitstat_total; merge fitstat_total fitstat_k&r; by Descr; run; data dim_k&r; set dim_k&r; rename Value=Value_K&r; run; /* Add to dataset dim_total with dimensions from each K=r */ data dim_total; merge dim_total dim_k&r; run; %end; /* End looping over values of K */ /****************/ /* Print Results of Interest */ /****************/ ods listing; data ae_total; set ae_total; rename Label=Parameter; run; data parm_total; set ae_total parm_total ; run; data fitstat_total; set fitstat_total; rename Descr=Criteria; run; /* Determine which K is selected by the preferred information criteria */ %if &info_crit="BIC" %then %do; data fit_crit; set fitstat_total; where Criteria="BIC (smaller is better)"; run; %end; %if &info_crit="AIC" %then %do; data fit_crit; set fitstat_total; where Criteria="AIC (smaller is better)"; run; %end; proc transpose data=fit_crit out=tfit_crit; by criteria; run; data tfit_crit; set tfit_crit; K_number=(_N_ -1); run; proc sort data=tfit_crit; by COL1; run; data tfit_crit; set tfit_crit; rank = _N_; run; data tfit_crit; set tfit_crit; where rank=1; run; proc sql noprint; select K_number into: K_select from tfit_crit; quit; %if &all_K=0 & &K_select=1 %then %do; %let K_select=&Kmax; %end; /* For the K selected combine the parameter estimates and additional estimates in one dataset */ %do r=&K_select %to &K_select; data ae_k&r; set ae_k&r; rename Label=Parameter; run; data parm_final_K; set ae_k&r parm_k&r; run; %end; title "Model Information"; title2 "By Degree of Flexibility"; proc print data=dim_total noobs; run; title "Information Criteria"; title2 "By Degree of Flexibility"; proc print data=fitstat_total noobs; run; title "Parameter Estimates"; title2 "By Degree of Flexibility"; proc print data=parm_total noobs; run; title "Parameter Estimates"; title2 "K = &K_select"; title3 "K selected by"; title4 &info_crit; proc print data=parm_final_K noobs; run; /*******************/ /* Graph estimated densities */ /*******************/ axis1 label=("Random Coefficient 1"); axis2 label=("Random Coefficient 2"); title "K=0 Estimated Random Effects Density"; proc gcontour data=graphics; plot x_value2*x_value1=density_K0 / haxis=axis1 vaxis=axis2 join pattern; run; %do r=(1*&all_K+&Kmax*(1-&all_K)) %to &Kmax; title "K=&r Estimated Random Effects Density"; proc gcontour data=graphics_k&r; plot x_value2*x_value1=density_k&r / haxis=axis1 vaxis=axis2 join pattern; run; %end; quit; %end; /*End q=2 Code*/ %mend;