9 Comments

1. 

   Dale McLerran on November 20, 2017 1:39 pm

   For calculating the pmf for displaying expected frequencies in the histogram, you assume that the number of Bernoulli trials AND the parameters a and b are known quantities. Of course, when presented with data believed to come from a beta-binomial distribution, the parameters a and b will typically need to be estimated. Even though specified here as a parameter of the beta-binomial distribution, the number of Bernoulli trials will generally be known. When the number of trials are known, the parameters a and b can be estimated. The Wikipedia entry on the beta-binomial distribution mentions both moment estimators of a and b and also states that one can employ likelihood maximization to estimate the unknown parameters. Note that estimation through likelihood maximization allows the parameters themselves to depend on predictor variables . We do not have predictor variables in the example provided. But should predictor variables exist, the ability to allow the parameters to depend on covariates is beneficial.

   Below, I present NLMIXED code that produces ML estimates of the parameters a and b. The NLMIXED procedure also computes the expected frequency under the beta-binomial distribution. There is no need for a data step to produce the expected frequencies. The parameter estimates are read in a data step and written to macro variables that are subsequently referenced as &a. and &b. in the SGPLOT procedure.

   Note that the default NLMIXED likelihood maximization procedure (QUANEW) produced estimates of the parameters but did not produce standard errors for the parameter estimates. For these data, the DBLDOG maximization technique produced estimates with confidence intervals that contain the true parameters. Other data may require a different likelihood maximization technique.

   Note, too, that the parameters a and b must be positive. Consequently, the NLMIXED procedure is parameterized in terms of log_a and log_b. Those parameters are exponentiated to produce a and b which are passed into the log-likelihood computation. Should a and b depend on covariates, the code expressed as "a = exp(log_a)" could be rewritten as "a = exp(b0_a + b1_a*X1 + b2_a*X2 ...;" with similar code for parameter b.

   ods output AdditionalEstimates=betaBin_ests; proc nlmixed data=BetaBin tech=dbldog; parms log_a log_b 1; /* Parameterize in terms of log_a and log_b */ a = exp(log_a); /* a > 0 ensured */ b = exp(log_b); /* b > 0 ensured */ nTrials = 10; ll = lcomb(nTrials, x) + logbeta(x+a, nTrials - x + b) - logbeta(a, b); pmf = exp(ll); model ll ~ general(ll); estimate "a" exp(log_a); estimate "b" exp(log_b); predict &SampleSize. * pmf out=PDFBetaBinom_est(keep=x pred); /* expected frequency computation */ run;   /* Keep just a single expected value record for each unique value of X. All other */ /* records are duplicates. Also, the predict statement names the predicted value */ /* as "pred". Rename to EX_est for consistency with Rick Wicklin presentation. */ proc sort data=PDFBetaBinom_est(rename=(pred=EX_est)) nodupkey; by x; run;   /* Write parameter estimates a and b to macro variables */ data _null_; set betaBin_ests; call symput(label, put(estimate,6.4)); run;   /* Merge simulated data and PMF. Overlay PMF on data distribution. */ data All_est; merge BetaBin PDFBetaBinom_est(rename=(x=t)); run;   title "Estimated Beta-Binomial Distribution"; title2 "Sample Size = &SampleSize"; proc sgplot data=All_est; vbarbasic x / barwidth=1 legendlabel='Simulated Sample'; /* requires SAS 9.4M3 */ scatter x=t y=EX_est / legendlabel='Expected Value' markerattrs=GraphDataDefault(symbol=CIRCLEFILLED size=10); inset "nTrials = 10" "a = &a." "b = &b." / position=topleft border; yaxis grid; xaxis label="x" integer type=linear; /* force TYPE=LINEAR */ run;

   Reply
   • 

     Rick Wicklin on November 20, 2017 2:14 pm

     Thanks for writing and for including the code for the MLE estimates. I like the trick of using the parameters on the log scale. Readers who want to learn more about how to use PROC NLMIXED for MLE can refer to the article "Two ways to compute maximum likelihood estimates in SAS."

     Since you mentioned regression, I will mention that PROC FMM can support beta-binomial regression. It uses a different parameterization and does not directly estimate a and b, but only the mean value a/(a+b), which is often the quantity of interest. It also requires that the data be in events/trials syntax. If you keep the nTrials variable in the simulated data, you can run the following:

     proc fmm data=BetaBin; model x/Ntrials = / dist=betabinomial; output out=out predicted; run;   data _null_; set out(obs=1); a = 6.8352; b = 4.5073 ; mu = a/(a+b); diff = Pred - mu; put Pred= mu= Diff=; run;

     Reply
     • 

       Travis on March 23, 2018 10:47 am

       Rick. I'm trying to find help, getting beta (PROC GLIMMIX Beta distribution) back on original scale

       PROC glimMIX;
       CLASS Phase TRT DAY ID;
       MODEL DMIBWnew = TRT|phase /dist=beta DDFM=KR SOLUTION;
       Random DAY/residual SUBJECT=ID;
       LSmeans trt /DIFF ADJUST=SIMULATE (REPORT SEED=121211) cl adjdfe=row;
       RUN; Quit;

       Reply
       • 

         Rick Wicklin on March 23, 2018 10:51 am

         I see that you have posted your question to the SAS Support Community. That's the correct place to ask questions like this.

         Reply
2. 

   Dale McLerran on November 20, 2017 6:08 pm

   The NLMIXED code can be rewritten to use a parameterization like that employed by FMM. See the code below. Of course, FMM would be the preferred procedure most of the time. But if one needed to introduce random effects due to data clustering, one could modify the NLMIXED code to accommodate that. Note that this parameterization appears more robust than the previous parameterization in that the default likelihood estimation technique (QUANEW) produced estimates with standard errors. (The DBLDOG technique also returns estimates with standard errors.)

   Note that parameters a and b in this formulation of the likelihood match parameters from the previous formulation except at the 5th significant digit of b. Moreover, the estimate for eta matches the intercept estimate from FMM through the 4 significant digits reported by each procedure. The scale (overdispersion) parameter differs by a single digit at the 6th significant digit.

   ods output AdditionalEstimates=betaBin_estsV2; proc nlmixed data=BetaBin; parms b0 0 log_a 1; eta = b0; pi = 1 / (1 + 1/exp(eta)); a = exp(log_a); phi = a/pi; nTrials = 10; mu = nTrials*pi; ll = lgamma(nTrials+1) - lgamma(x+1) - lgamma(nTrials - x + 1) + lgamma(phi) - lgamma(nTrials+phi) + lgamma(x+a) + lgamma(nTrials - x + phi - a) - lgamma(a) - lgamma(phi-a); model ll ~ general(ll); estimate "a" a; b = a*(1-pi)/pi; estimate "b" b; estimate "scale" phi; pmf = exp(ll); predict &SampleSize. * pmf out=PDFBetaBinom_est(keep=x pred); /* expected frequency computation */ run;

   Again, since FMM can estimate this same beta-binomial model, the only reason to keep the NLMIXED code in one's hip pocket would be if there were yet another source of overdispersion attributable to random effects associated with clustered data collection.

   Reply
   • 

     Rick Wicklin on November 21, 2017 5:18 am

     Very nice summary. Thanks for taking the time to share this information.

     Reply
3. 

   Ksharp on November 25, 2017 7:02 am

   Rick and Dale,
   How could I apply this Beta-Binomial distribution into Logistic Regression ?

   Reply
4. 

   KAUSHIK BHATTACHARJEE on December 23, 2017 2:51 am

   Hi Rick,

   suppose Yt=f(Xt ,Zt) .Here only Yt is observable but Xt and Zt are not observable-however both are random variables for a given t( and hence{Xt} and {Zt} can be considered as two unobservable time series .)Therefore the problem is to disentangle the effect of Xt and Zt from Yt

   like KALCVS or KALDFF any subroutine or proc like SSM forPARTICLE FILTERING?

   Best,
   Kaushik

   Reply
   • 

     Rick Wicklin on December 23, 2017 7:12 am

     You can ask questions like this on the SAS Support COmmunities. Since this is related to time series, I suggest you post to the Forecasting and Econometrics Community.

     Reply