I have a macro for performing HLM and logistic regression and then outputting a results table with the odds ratios and 95% confidence intervals from both models. When I apply the macro and try to run it, it is not able to generate the upper 95% CI for the logistic regression (OR_lr95UCI). If I remove this line from the code, it generates a table with all other results but when I include this variable, the code does not run. I am calculating it the same way as the lower limit (OR_lr95LCI), so I don't know what is causing the issue. I have copied the macro below, with the problem variable bolded. %MACRO HLM2(data=, outcome=, covariates=, data2=, covariates2=, b2index=); /* Run PROC LOGISTIC, result is in the data set logistic_output */ proc logistic data=&DATA desc COVOUT OUTEST=logistic_output; model &OUTCOME = &COVARIATES ; run; proc iml; /**************************************************************** *** Hirearchical model using empirical-Bayes or semi-Bayes estimation. *** ****************************************************************/ start hm_reg(b,v,z,t2,max_iter,c_criter,b_hm,bcov_hm,b2); /**************************************************************** ** Input: b : 1st-stage parameter estimates ** v : 1st-stage covariance matrix estimates for b ** z : 2nd-stage design matrix (prior design matrix) ** t2 : prior variance for semi-Bayes estimation (vector), ** (0 if prior variance is to be estimated using ** emirical-Bayes estimation.) ** max_iter: maximum number of iteration (default is 20) ** c_criter: convergence criteria (default is 1e-3) ** Output: b_hm : posterier estimates ** bcov_hm : covariance matrix for b ** b2 : 2nd-stage parameter estimates ****************************************************************/ /*** Number of first-stage parameters.***/ np = nrow(Z); ncolz = ncol(Z); /*** Second-stage degrees of freedom. */ df2 = np-ncolz; Inp = I(np); /*** Empirical-Bayes estimation, else semi-Bayes estimation. ***/ if (min(t2) <= 0) then t2pass = 0; else t2pass = 1; /*** Maximum number of iterations, and convergence criterion for t2. ***/ if max_iter = 0 then _maxiter = 20; else _maxiter = max_iter; if c_criter = 0 then _criter = 1e-3; else _criter = c_criter; if t2pass = 0 then do; print np[format=4.0] " first stage parameters were regressed on" ncolz[format=4.0] " second stage covariate(s)", " using empirical-Bayes estimation method (prior variance (t2) to be estimated)."; end; else do; print np[format=4.0] " first stage parameters were regressed on" ncolz[format=4.0] " second stage covariate(s)", " using semi-Bayes estimation with prior variance (t2) fixed at" t2[format=8.4]; end; /*** Undertake second-stage regression: ***/ t2old = t2+1; do iter = 1 to _maxiter while(abs(t2-t2old) > _criter); t2old = t2; /*** Weight matrix. ***/ w = inv(v+t2#Inp); wv = w*v; wz = w*z; ws = sum(w); /*** Invert 2nd-stage information. ***/ vs = inv(t(Z)*w*Z); /*** 2nd-stage coefficient estimates. ***/ b2 = vs*(t(wz)*b); if t2pass=0 then do; /*** Residual from 2nd-stage estimates. ***/ e = b - (Z*b2); /*** Total residual sums of squares. ***/ rsst = t(e)*w*e; /*** Residual mean square. ***/ rms = np*rsst/(df2*ws); /*** Update t2 for empirical Bayes. ***/ t2 = max(rms-sum(wv)/ws,0); end; end; /*** Calculate posterior expectations of 1st-stage parameters: ***/ /*** Weight for prior expectations of the 1st-stage parameters (curvature correction for EB.) ***/ if t2pass=0 then wvc = ((df2-2)/df2)*wv; else wvc = wv; /*** Projection of b to prior mean. ***/ hatw = Z*vs*t(wz); st2 = sqrt(t2*t(t2)); /*** Projection of b to posterior mean. ***/ hatp = wvc*hatw + w#st2 + (wv-wvc); /*** Estimated posterior mean. ***/ b_hm = hatp*b; /*** Estimated posterior covariance. ***/ bcov_hm = v - wvc*(Inp-hatw)*v; /*** Approx df in hierarchical model. ***/ npa = trace(hatp); if t2pass=0 then do; /*** Add variance component due to estimation of the prior variance: ***/ wvce = wvc*(e#(sqrt((sum(wv)/ws + t2)/(vecdiag(v)+t2)))); bcov_hm = bcov_hm + (2/(df2-2))*wvce*t(wvce); /*** Check for fitting problems: ***/ if iter >= _maxiter then print " Empirical-Bayes did not converge."; else do; if t2=0 then print " Empirical-Bayes converged after" iter[format=4.0] " iterations," " but potential underdispersion: t2 set to zero."; else print " Empirical-Bayes converged after" iter[format=4.0] " iterations," " with prior variance (t2) estimate:" t2[format=8.4]; end; end; do i=1 to np; if bcov_hm[i,i]<0 then do; print "Negative second-stage variance: set to zero."; bcov_hm[i] = 0; end; end; finish; /****************************************************************/ /**************************************************************** *** Main procedure for the estimation. Check model setting for each run *** ****************************************************************/ /************************* **** Ordinary Logistic Regression Model *** *************************/ bnames = {INTERCEPT &COVARIATES}`; /* Make use of logistic_output */ use logistic_output; read all into lr_est; b_lr = lr_est[1,1:(ncol(lr_est)-1)]`; bcov_lr = lr_est[2:nrow(lr_est),1:(ncol(lr_est)-1)]; /********************************* *** Output logistic result *** *********************************/ SE_lr = sqrt(vecdiag(bcov_lr)); ward = (b_lr # b_lr) / (se_lr # se_lr); P_value = 1 - probchi(ward,1); print "Ordinary logistic regression estimates:"; print bnames b_lr[format=10.4] se_lr[format=10.4] p_value[format=8.4]; /******************* *** HLM Model *** *******************/ /*** Initialize 1st-stage parameters to be modelled in the hierarchical model ***/ b_lr0 = b_lr[&b2index]; v_lr0 = bcov_lr[&b2index,&b2index]; bnames_hm = bnames[&b2index]; /*** Specify 2nd-stage data and parameters to be modelled in the hierarchical model ***/ use &DATA2 var {&COVARIATES2}; read all into Z_mat; bnames2 = {&COVARIATES2}`; print bnames2 Z_mat; Z_mat = repeat(1,nrow(b_lr0),1); bnames2 = {INTERCEPT}; t2 = 0.5; run hm_reg(b_lr0, v_lr0, Z_mat, t2, 0, 0, b_hm, bcov_hm, b_hm2); /**************************** *** Output HLM result *** ****************************/ SE_hm = sqrt(vecdiag(bcov_hm)); ward = (b_hm # b_hm) / (SE_hm # SE_hm); P_value = 1 - probchi(ward,1); print "Hierarchical logistic regression estimates:"; print bnames_hm b_hm[format=10.4] SE_hm[format=10.4] P_value[format=8.4]; /*** OR estimates from both ordinary logistic model and hierarchical model ***/ OR_lr = exp(b_lr[&b2index]); OR_lr95LCI = exp(b_lr[&b2index] - 1.96*SE_lr[&b2index]); OR_lr95UCI = exp(b_lr[&b2index] + 1.96*SE_lr[&b2index]); OR_hm = exp(b_hm); OR_hm95LCI = exp(b_hm - 1.96*SE_hm); OR_hm95UCI = exp(b_hm + 1.96*SE_hm); /*** Print results: ***/ print "Second stage parameter estimates using hierarchical model.", bnames2 b_hm2[format=10.4]; print " Odds Ratios and 95% CIs", " Conventional Model Hierarchical Model", bnames_hm OR_lr[format=8.2] OR_lr95LCI[format=8.2] OR_lr95UCI[format=8.2] OR_hm[format=8.2] OR_hm95LCI[format=8.2] OR_hm95UCI[format=8.2]; quit; run; %MEND;
... View more
