I have the following simulation code to produce 1000 sets of data similar to "work.final" where the independent variables do not change in all the 1000 sets.   But I would like to add another condition. For each of the 1000 datasets, the frequency of the value 10, which is a simulated value of y should be 20%, the same way...the frequency of value 15 should be 10%.   %let NumSamples = 1000; data s_new; call streaminit(6578); set work.final; do SampleID=1 to &NumSamples; mu = exp(1.5360 + 0.0140*x1 - 0.1016*x2 + 0.6216*x3); k = 1/0.5345; p = 1/(1+0.5345*mu); y = rand('negbinomial', p, k); output; end; run;   ... View more

I have written a procedure that uses a specific estimate to compute Poisson probabilities for values between 0 to 80. But at specific values (5, 12, 20, 30, 45, 50, 60), the procedure for computing the probability is a little different as in the code.    DATA xxxx2;   SET final_sim_data;   estimate=3.99793;   IF x1~=1 or x2~=1 or x3~=1 then DELETE; run; data pred; set xxxx2; array v[11]; array n[11]; array sf[11]; array t[11]; array f[11]; array a[11]; array w[11]; do i = 1 to 11; n[i] = i + 14; end; do i = 1 to 11; sf[i] = i + 44; end; do i = 1 to 11; v[i] = i + 54; end; do i = 1 to 11; t[i] = i + 24; end; do i = 1 to 11; f[i] = i + 39; end; do i = 1 to 11; a[i] = i - 1; end; do i = 1 to 11; w[i] = i + 6; end; do count_y= 0 to 80; if count_y ~= 5 and count_y ~=12 and count_y ~=20 and count_y ~=30 and count_y ~=45 and count_y ~=50 and count_y ~=60 then probability_FP= pdf("poisson",count_y,exp(estimate)); else if count_y= 20 then do; p = 0; do i = 1 to 6; p = p + ((n[i]-15)/5)*(pdf("poisson",n[i],exp(estimate))); end; do i = 7 to 11; p = p + ((25-n[i])/5)*(pdf("poisson",n[i],exp(estimate))); end; probability_FP=p; end; else if count_y= 60 then do; sum = 0; do i = 1 to 6; sum = sum + ((v[i]-55)/5)*(pdf("poisson",v[i],exp(estimate))); end; do i = 7 to 11; sum = sum + ((65-v[i])/5)*(pdf("poisson",v[i],exp(estimate))); end; probability_FP=sum; end; else if count_y= 30 then do; th = 0; do i = 1 to 6; th = th + ((t[i]-25)/5)*(pdf("poisson",t[i],exp(estimate))); end; do i = 7 to 11; th = th + ((35-t[i])/5)*(pdf("poisson",t[i],exp(estimate))); end; probability_FP=th; end; else if count_y= 50 then do; sd = 0; do i = 1 to 6; sd = sd + ((sf[i]-45)/5)*(pdf("poisson",sf[i],exp(estimate))); end; do i = 7 to 11; sd = sd + ((55-sf[i])/5)*(pdf("poisson",sf[i],exp(estimate))); end; probability_FP=sd; end; else if count_y= 12 then do; bb = 0; do i = 1 to 6; bb = bb + ((w[i]-7)/5)*(pdf("poisson",w[i],exp(estimate))); end; do i = 7 to 11; bb = bb + ((17-w[i])/5)*(pdf("poisson",w[i],exp(estimate))); end; probability_FP=bb; end; else if count_y= 5 then do; aa = 0; do i = 1 to 6; aa = aa + ((a[i]-0)/5)*(pdf("poisson",a[i],exp(estimate))); end; do i = 7 to 11; aa = aa + ((10-a[i])/5)*(pdf("poisson",a[i],exp(estimate))); end; probability_FP=aa; end; else if count_y= 45 then do; fr = 0; do i = 1 to 6; fr = fr + ((f[i]-40)/5)*(pdf("poisson",f[i],exp(estimate))); end; do i = 7 to 11; fr = fr + ((50-f[i])/5)*(pdf("poisson",f[i],exp(estimate))); end; probability_FP=fr; end; output; end; keep count_y probability_FP; run; /********************/ proc means data=pred; class count_y; output out=means mean(probability_FP)=; run;     1. I am sure the code is so manual (I mean it could be shortened) especially because I am averaging the same value in the last proc code. Is there a way I can shorten everything so that I avoid the proc means code?   2. After getting the probabilities for count_y = 0,1,2...80. I notice that the probabilities are not summing to 1. What could be the problem, and Is there a way out. ... View more

I want to simulate Zero inflated Poisson data sets 1000 times. I am using predicted coefficients and independent variables that are set as in the code below. But I have to include more conditions;   1. I want 70% of the resulting counts(the simulated y) in the interval [5,14] be rounded of to 10 for each data set.       The thing is that the first code below simulates each observation 1000 times, then i'm using the second code to order the data so   that I have 1000 data sets. So the rounding is required for each data set.   %let NumSamples = 1000; /* number of samples */ data Sim_data(drop=mu); set WORK.xxx; call streaminit(1234); do SampleID=1 to &NumSamples; ObsNum = _N_; mu = exp(0.6773 - 0.01036*Age + 0.2080*height); ypoi = ranpoi(1234,mu); pzero = cdf("LOGISTIC",-0.9792 +0.03796*Age -0.1488*height); if ranuni(1234)>pzero then do; ypoizim = ypoi; end; else do; ypoizim = 0; end; y=ypoizim; output; end; run;       proc sort data=Sim_data; by SampleID ObsNum; run; 2. Again the same thing, rounding of 20% of resulting counts in the interval [25,35] to 30 ... View more

  I'm running the following two codes on the attached data and getting two different results. Moreover the results for the Zero inflation do not appear correct. I don't know whether it is a problem with some code or It's a sign that zero inflation is not the best for this data! because even the significance of the Zero inflation estimates is way off from proc nlmixed     proc genmod data = bbbb; model y_c = x1 x2/dist=zip; zeromodel x1 x2 /link = logit; run; /*********************************************************************************************************/ proc nlmixed data=bbbb; parameters b0=0 b1=0 b2=0 a0=0 a1 = 0 a2=0 ; linpinfl = a0 + a1*x1 + a2*x2; infprob = 1/(1+exp(-linpinfl)); lambda = exp(b0 + b1*x1 + b2*x2 ); if y_c=0 then ll = log(infprob + (1-infprob)*exp(-lambda)); else ll = log((1-infprob)) - lambda + y_c*log(lambda) - lgamma(y_c + 1); model y_c ~ general(ll); run; ... View more

I have two models estimating the same number of parameters  and using the same data. The models are only different because they are using different different likelihood functions in the proc nlmixed procedure. I mean this;   proc nlmixed data=Surveydata; parms b0=0 b1=0 ; mean = exp(b0 + b1*Age); ll = /*Some codes*/; model art ~ general(ll); run; proc nlmixed data=Surveydata; parms b0=0 b1=0 ; mean = exp(b0 + b1*Age); ll = /*Some codes different from those above*/; model art ~ general(ll); run; Therefore expect different estimates for variable Age. So my question is this Apart from AIC or BIC what test can I use especially that I can incorporate with simulation to justify that one model is better than the other. The test may even focus on the ability of the estimates of the parameter Age to best explain the variable art. But most of all I would really want to simulate different data sets to use the test just to justify that for different data sets Model X is better than Model 2.   ... View more