When you coded the log likelihood, this might result if you used LOG(GAMMA(...)) instead of the LGAMMA function as shown in the examples I referred to earlier. ... View more

See this note on overdispersion. Your Poisson model suggests overdispersion as evidenced by the Pearson or Deviance divided by the DF being considerably larger than 2. The negative binomial is often used for overdispersed data. If you do that, these statistics are only about 1.1 suggesting that overdispersion is not a problem with that model. Or, as suggested in the note, another common approach is to use a GEE model which seems to give a similar estimate on your primary predictor as the non-GEE negative binomial model (using either Poisson or negative binomial distribution in the GEE model). proc genmod data=forSAS namelen=65; class Female Location Income_thresh id; model Errors=LPA_Daily_min ST_win Female Location Income_thresh / dist=negbin; repeated subject=id; effectplot fit(x=LPA_Daily_min) / ilink; run; ... View more

The log likelihood for the beta-binomial is shown in "Log-Likelihood Functions for Response Distributions" in the Details section of the FMM documentation. As suggested, you can use PROC NLMIXED to fit a model specifying this log likelihood and adding a RANDOM statement if desired. You can see examples of fitting a model with specified log likelihood function in this note showing this using the log likelihood for the truncated Poisson and negative binomial distributions and in this note where the zero-inflated Poisson and negative binomial are specified. ... View more

You can easily get a plot showing the effect of a continuous predictor on the event probability using an EFFECTPLOT statement. For example, adding this statement in the Getting Started example in the PROC GEE documentation gives a plot of the fitted model showing the effect of Age with other variables fixed. effectplot fit(x=age)/ilink; ... View more

The lack of balance, meaning unequal numbers of observations in the various predictor combinations, is not itself a problem. It is the extreme case of this when some of the combinations have no observations. That is when "separation" occurs resulting in some parameters being infinite as I mentioned. In these cases, you might be able to fit the model with at least some interactions by using a penalized likelihood. This can be done by simply adding the FIRTH option in the MODEL statement. Regarding your question - you don't need two WHERE statements because you can specify a single WHERE statement with multiple conditions such as: where varietey='Hass" and Season=2022; ... View more

Whether you use a mixed model with GLIMMIX or a Generalized Estimating Equations model with PROC GEE, you can get odds ratio estimates for your categorical (in  CLASS statement) by specifying them in an LSMEANS statement with the DIFF and ODDSRATIO options. If you specifically want to estimate relative risks, then see the available methods discussed in this note. ... View more

Actually, I don't think this will help - the degrees of freedom of the other interactions are absorbed into the 3-way so you end up with just as many parameters requiring estimation. ... View more

The problem is not the zero values in the response variable. For a multinomial categorical variable, zero is just another category and the distribution does not restrict the proportion of zeros like with, say, the continuous gamma distribution. The problem here is that you have several response categories and also specify all possible interactions resulting in a complex model with many parameters to be estimated and the model complexity makes the data in each variety too sparse. The result, just like for binary logistic models, is that some model parameters are actually infinite. Since computers don't deal in infinities, the practical result is that some parameters are large with even larger standard errors and/or some parameters have zero degrees of freedom and are not estimated. The solution is to simplify the model in any way acceptable to you such as by combining response levels (which will have the biggest benefit) and/or removing some or all interactions. Because the amount of data varies so much by variety, you should not expect that you can just specify one model and use BY VARIETY and get a proper fit for each variety unless you find a much simpler model that can be successfully fit in every variety. And by the way, for logistic models like this (binary or multinomial), PROC LOGISTIC is the better procedure to use than PROC GENMOD as it is more specialized for those models. ... View more

It uses the last level when internally sorted in alphanumeric order. So, with levels A, B, and C, it will use C as the reference level. ... View more

Yes. Hundreds of clusters would generally be considered adequate for validity of the method. ... View more

There are several types of models that can be used to model clusters of correlated observations like you describe. One easy approach is the Generalized Estimating Equations (GEE) model which you can fit using PROC GEE. See the example for modeling a binary response (like your death response) using this model in the Getting Started section of the PROC GEE documentation. This model does not require equal numbers of observations per cluster (subject) nor that observations be at the same times, particularly if you use a correlation structure like TYPE=IND or EXCH as in that example. The predictors can vary over time within subjects. You should, however, have a large number of clusters for validity of the method. ... View more

But since the risk ratios you want are just comparisons among particular pairs of levels, the equivalent can be done without the ordinal parameterization using the easier LSMEANS method shown in the example in the note.  proc sort data=neur; by descending duration; run; proc genmod data=Neur; class Treatment Sex id duration(order=data); model painnum= Treatment Sex Age Duration / d=p; repeated subject=id; lsmeans duration / diff exp cl; run; ... View more

The first link I provided in my post shows the Poisson method using PROC GENMOD that you're talking about. In the example there, the CLASS predictor of interest is treated as nominal and in that case, the easiest way to get risk ratios is using GLM parameterization and the LSMEANS statement as shown rather than to mess with ESTIMATE statements.    In your case using the ORDINAL parameterization, the parameters of the variable are interpreted as the incremental change resulting from increasing the variable from the previous level to each level in turn.  So, the relative risk associated with each level change is just the parameter exponentiated. Using the neuralgia data in the example titled "Logistic Modeling with Categorical Predictors" in the PROC LOGISTIC documentation, these statements create an ordinal CLASS variable, DURATION, fit the Poisson model with it using ORDINAL parameterization, and use ESTIMATE statements to produce the relative risks which appear in the Mean columns.  proc rank data=neuralgia groups=4 out=neur; var duration; run; data neur; set neur; id=_n_; painnum=(pain="Yes"); run; proc genmod data=neur; class Treatment Sex duration(param=ordinal) id / param=ref; model painnum= Treatment Sex Age Duration / d=p; repeated subject=id; estimate 'rr 0-1' duration 1 0 0; estimate 'rr 1-2' duration 0 1 0; estimate 'rr 2-3' duration 0 0 1; run; ... View more

So you do apparently want statistical power. Power and sample size analysis can be done in PROC POWER for a planned study to be analyzed using a logistic model. See the example in the PROC POWER documentation. Note that PROC POWER is for prospective power analysis, not for a retrospective study that has already been designed and analyzed. Results from such a prospective analysis done in advance of you collecting data is apparently what your reviewer is asking for. ... View more