Check that note out -that an R side variance component is confounded with the error.   So go ahead and identify a row-level ID, which hopefully is applicable across observational IDs. Call it rowblah, I guess. Then try the following (no guarantees that the whole thing won't blow up completely):   proc glimmx data=input empirical=mbn; class ID blah bla rowblah; logtime=log(pyrs); effect aspl=spline(cov1 / NATURALCUBIC BASIS=TPF(NOINT) notmethod=PERCENTILELIST(5 27.5 50 72.5 95)); effect pspl=spline(cov2/ NATURALCUBIC BASIS=TPF(NOINT) knotmethod=PERCENTILELIST(5 27.5 50 72.5 95)); *class ID blah someproperty(ref='0'); model events=blah aspl pspl bla rowblah aspl*rowblah pspl*rowblah/ dist=poisson offset=logtime s cl; random rowblah / residual subject=ID; run; I commented out the second CLASS statement, as I think it is redundant, and someproperty does not appear in either the MODEL or RANDOM statement.  This approach essentially treats each level of rowblah as a repeated measure on the subject.  I don't know if it will work, as convergence issues, starting value issues, and other possible problems may be lurking.  However, it strikes me as the most easily implemented method.   SteveDenham   ... View more

I am going to assume that you are using the brms package in R.  Good but takes a long time, especially with as many levels as you have.  I suspect the start-up/burn-in stage suffers from the issue you are seeing with GLIMMIX.  However, the Bayes approach uses the starting values supplied and generates new realizations at each iteration.  That makes me think that maybe, just maybe, using PROC MCMC in SAS might work a bit faster.  See if Fang Chen's SGF paper here https://support.sas.com/resources/papers/proceedings16/SAS5601-2016.pdf  for more on this.   SteveDenham ... View more

I vote Bayesian. Didn’t know that was an option. I’ll ask some things about it tomorrow. SteveDenham ... View more

I really sympathize with your desire to "show up" the others in your field that consolidated low count categories.  That would be great.  However, the first thing that occurs to me is that they may have run into the same issue, and took that alternative as the only way to get an analytic solution.   I guess the only thing I could suggest at this point is to see if there is any way to turn this into an ordinal multinomial. Unfortunately I don't have enough knowledge in this field to even begin to suggest how that might be done.  Tell us something about the independent variable TrgPos.  Could the categories there be assigned some sort of ordinal value, based on some utility function?  If so, then you could move from the glogit link and many of the issues you have seen to a clogit link.  The random effects wouldn't change, but at least it would be something different.   SteveDenham ... View more

It is the problem, in the sense that you present it.  I suppose if you had terabytes of CPU that wouldn't be a problem.   What can be done? First, you need to consolidate some of these.  You have region, family within region, and individual within family (=residual).  An R side approach in GLIMMIX may be useful (but I worry about twins under this approach). Consider family as a repeated measure within region.  If you aggregate at the family level, and assign a weight = number of family members, this might work (no guarantees, though):   Proc glimmix data=temp; class toxin yeargrp SES f_region ; nloptions maxiter=500 tech=nrridg; model case=toxin yeargrp SES f_region/ dist=poisson link=log offset=lnpyrs covb cl solution; random f_region/residual subject=family_id type=cs; weight = famnumbers; run; There are some important considerations here.  First, you need a unique numeric family_id for each family so that it can be treated as a continuous variable,  It would be a good idea to sort the dataset by this variable. Second, toxin yeargrp and SES may have different effects in each region.  This model would give the marginal effects averaged over region.  To get effects for each region, you would have to specify interactions. Third, case will have to be aggregated on a family level.  This part is not too difficult using PROC MEANS, which would also give the number of observations for each family (famnumbers).  The join on the design matrix (original dataset) may be difficult.   My fear is that if the clustered SD approach failed in GENMOD due to memory constraints, then it will happen here as well.  Jack-knifing the big dataset into smaller sets that will run, and then consolidating the results, might be an approach to consider then.   SteveDenham ... View more

No you are not.  Just as I replied to @PaigeMiller regarding dep_var in the CLASS statement, this is another quirk of nominal multinomial analysis in GLIMMIX.  I don't see anything about it in the documentation (could be looking in the wrong place), but it surely makes sense. My suggestion was not useful, and I apologize.   So, now we come to the issue of the random variables CompLemma and trigger.  What can you say about these?  How many levels are there in each?  I note that they are not nested in any manner, but that is about all I know.  You said that if CompLemma is the only random effect specified that you get a message that the G matrix is not positive definite.  That hints at something within that variable being wonky.  Try your original code with a NOBOUND option in the PROC GLIMMIX statement.  It may be one of those situations that in an OLS setting would lead to a "negative variance component".  That means it is strongly correlated with trigger (negatively) or the residual error, which in turn could mean you have a level or levels that could be consolidated.  Look at some PROC FREQ output for these variables.  If you have levels with small N (and that is going to be subjective), consider consolidating those levels, provided that makes sense for your research questions.  We just need more information about these variables, and why they need to be considered.   One final approach might be to run a model with CompLemma and trigger as fixed effects.  Perhaps the solution vector would give some estimates that could be used as starting values in a PARMS statement.  I don't know???   SteveDenham ... View more

When you say "does not converge", the first thing I look for is an NLOPTIONS statement, where you can increase the number of iterations beyond the default 20.  For now try adding     nloptions maxiter=1000; If you still do not converge, then look at this paper for some really good ideas about improving convergence in mixed models:   https://support.sas.com/resources/papers/proceedings12/332-2012.pdf    SteveDenham   ... View more

That should not present a problem - once again it is why LSMESTIMATE statements are so useful.   lsmeans vaper|riagender/e; lsmestimate riagender "Males vs Females, averaged over vaper status' -1 1/exp; lsmestimate vaper*riagender 'Males vs Females, vaper status=never' -1 1 0 0 0 0, 'Males vs Females, vaper status=former' 0 0 -1 1 0 0, 'Males vs Females, vaper status=current' 0 0 0 0 -1 1/exp; This now gets trickier.  I just stuck these in, and the correct version depends on the ordering of vaper levels.  Since there is nothing in the CLASS statement and I assume that the levels are in fact  'Current', 'Former', and 'Never', the lsmeans will sort in alphabetical order.  That means the code should look like:   lsmeans vaper|riagender/e; lsmestimate riagender "Males vs Females, averaged over vaper status' -1 1/exp; lsmestimate vaper*riagender 'Males vs Females, vaper status=current' -1 1 0 0 0 0, 'Males vs Females, vaper status=former' 0 0 -1 1 0 0, 'Males vs Females, vaper status=never' 0 0 0 0 -1 1/exp; If there is another coding for the vaper levels, it can be handled once you know the order that things appear in the LSMEANS statement.   SteveDenham.     ... View more

0.333 is the point estimate of the probability (2 out of 6).  The 95% exact confidence bounds are 0.0433 and 0.7772, so the true probability is likely to be in that interval in 95% of the calculations of the point estimate.   You could plug these bounds into the formulas I had earlier, and get estimates for the 50% fiducial limits (point limits) for the 2.5 %ile point and 97.5%ile point.   What i really think you might want is: Given a sample size of 60, what is the number of adverse events you would have to observe to be 95% confident that the true value was 0.333 with 80% power.  Now that, or something like that, is something you can plug into PROC POWER, using the onesamplefreq option.   SteveDenham ... View more

There is probably a way to do this in PROC POWER, but this is one case where I think a simple calculation can give the answer.  It is a matter of using the binomial formula [ NtakeM * pprob^M * (1-pprob)^(N - M) ].  So if you want the probability of at least one out of the 60, you calculate this as = (1 - prob(M=0)).  Using your prior estimate of pprob=0.333, I get prob(M=0) = 60take0 * (0.333)^0 * (1 - 0.333)^60 = 1 * 1 * 0.667^60 = 2.72e-11, which gives the probability of at least one adverse event as 1 - 2.72e-11 = 0.99999999997280278361063568173433.  The probability of exactly one adverse event is 60take1 * (0.333) * (1 - 0.333)^59 = 60 * 0.333 * 4.08e-11 = 8.1591649168092954797017275793104e-10.   SteveDenham ... View more

@Tidemann2nd , you nailed it.  Your code is the difference between the two group estimates.   SteveDenham ... View more

We do this on a regular basis using the LSMESTIMATE statement. If the variable of interest in the model is called treatment, then the following code would be usable:   lsmestimate treatment 'Early application' 1 0 1 0 1 0, 'Late application ' 0 1 0 1 0 1/divisor=3; This assumes equal weighting for each of the groups. The issue with using the format approach is that you would need a separate model with the formatted data to get the least squares means (I think).   SteveDenham   ... View more