Hi Rick, Thank you so much for your reply. You are right about the difference between prediction and simulation. I really want to simulate values based on the predicted (estimated) parameters. I did some checking on the simulated 0-1 values and it seems to be working fine. Simlar thing I want to do for normal districution also. That is to say, I would like to simulate from normal distribution (in two levels, as I am considering a mixed effects model) with different mean parameters but consatnat sd. If you please have a look at the folowing code (part 2) to assure that the code does exactly what I want to do, that would be very helpful. The code appears to be long but it should be relatively easier for you to go through it. /*Part1: Fit the logistic model to simulate synthetic 0-1 values*/ proc logistic data=nonfed noprint; class size naics2; model E_pos(event='1') = size naics2; output out=E_binary_pred p=phat_E; run; data Sim1; set E_binary_pred; E_pred_pos = rand("Bernoulli", phat_E); run; /*Part2: Fit the Linear Mixed Model for employment, given that it's positive*/ data E_positive; set nonfed; if E_pos = 0 then lAvgE=.; run; proc mixed data=E_positive method=reml; class size own naics2 state; id CvID EvID county naics; /*log employment(average of M1, M2, M3) has been modeled */ model lAvgE = size own naics2 MSAid / s OUTPM=E_LMM_XBetaHat;/*store the synthetic predicted values Xbeta_hat*/ random state / s subject=state type=VC; ods output CovParms = covparms_mixed_E;/*covparms_mixed_E contains variance component estimates*/ run; proc iml; use covparms_mixed_E; read all var {Estimate} into sigma2Est; close covparms_mixed_E; sigmav_hat = sqrt(sigma2Est[1]); sigmae_hat = sqrt(sigma2Est[2]); use E_LMM_XBetaHat; read all var {Pred} into mu; close E_LMM_XBetaHat; n = nrow(mu); sigmae = repeat(sigmae_hat, n); sigmav = repeat(sigmav_hat, n); use nonfed; read all var {AvgE AvgWage wages MSAid}; read all var {CvID EvID own size naics state county naics2}; close nonfed; create E_LMM_muHat var {CvID EvID own size naics state county naics2 AvgE AvgWage wages mu sigmav sigmae}; append; quit; data Sim2; set E_LMM_muHat; x2 = rand("Normal", mu, sigmav); x3 = rand("Normal", x2, sigmae); E_pred = exp(x3); run; data synth_E_selstate_datastep1; merge Sim1 Sim2; by EvID; run; data synth_E_selstate_datastep; set synth_E_selstate_datastep1; E_synth = E_pred_pos*E_pred; E_synth = round(E_synth,1); run; Thanks again, Santanu
... View more