I don't know which will give what you want, but I think that if you add 'female' and 'black' to your CLASS statement, the results should automatically handle the categorical nature of the variables.  It may be that this will give separate curves for each row in the 2x2 design (female x black), with the curves adjusted to the mean value within each row.  It would probably require fitting the interaction term to get separate curves for each cell.   But I am just sort of spitballing on this one.   SteveDenham ... View more

I can't guarantee that this is the issue but the syntax should be:   value = pdf('BINOMIAL', number_of_successes, probability_of_success, number_of_independent_trials)   Just looking at the values in the NOTEs, I would guess that you have switched the first and last values, as the line of data always has Number strictly greater than Ntarget.  Try exchanging the two variables in your call to the PDF function.:   SteveDenham     ... View more

A quick google on "combining log likelihoods" gives a lot of interesting results ranging from as simple as adding the log likelihoods (equivalent of multiplying the likelihood) which assumes independence between the imputations, to several Bayesian methods, and on to some good Q and A on stackexchange.   The one thing I thought was a common thread is that there is good theory behind the adding of log likelihoods/multiplying likelihoods when you think of it as combining probabilities.  There is a strong independence assumption though.   SteveDenham ... View more

NOTE: Convergence criteria met but final Hessian is not positive definite. NOTE: Estimated G matrix is not positive definite.   This is telling you that you have over-specified the random effects based on the data you have.   WARNING: Stopped because of too many likelihood evaluations.   This can be avoided in many cases by simply increasing the number of iterations using the MAXITER= option in the PROC MIXED statement.  However, look through the iteration history and see if it is approaching convergence, or if it is jumping around.  MAXITER will solve the first case, but the second will require more appropriate MODEL and RANDOM statements.   SteveDenham ... View more

The PROC POWER documentation is an excellent place to start for power and sample size calculations.  For distributions, it is hard to beat Merran Evans Statistical Distributions.   SteveDenham ... View more

I recommend looking into using LSMESTIMATE statements rather than ESTIMATE statements - the syntax is much more compact and intuitive.   SteveDenham ... View more

I like most of option 2.  I would consider fitting only the interaction term (with a NOINT option) and then using LSMESTIMATE statements (with an ILINK option) to get at the questions of interest.  Both option 2 and this method end up calculating marginal estimates for these effects - the proportions if all cells had an equal number of observations.  If you wish to get at the differences on the probability scale, then you will need the %NLmeans macro (and you should read the SAS Note(s) on this and everything @StatDave  has posted on this forum/topic.   SteveDenham ... View more

Try this:   proc glimmix data=WORK.IMPORT; Class NIVEL ENZ; Model CMMS = NIVEL ENZ NIVEL*ENZ/ddfm=kr; Lsmeans NIVEL*ENZ/ adjust=tukey diff lines e; contrast 'Linear NIVEL at ENZ=0' NIVEL*ENZ -3 -1 1 3 0 0 0 0; contrast 'Quadratic NIVEL at ENZ=0' NIVEL*ENZ 1 -1 -1 1 0 0 0 0; contrast 'Linear NIVEL at ENZ=1' NIVEL*ENZ 0 0 0 0 -3 -1 1 3; contrast 'Quadratic NIVEL at ENZ=0' NIVEL*ENZ 0 0 0 0 1 -1 -1 1; run; Since you are using CONTRAST rather than ESTIMATE, you can use integer coefficients.  If you want some output value, rather than just a p value, switch to LSMESTIMATE statements.   I added the 'E' option to your LSMEANS statement as a check on the arrangement of coefficients in the CONTRAST statements.   SteveDenham   ... View more

One thing to consider is the distribution of your response variable - continuous, count, proportion - as all need to handled a bit differently in PROC POWER.   SteveDenham ... View more

In point 2, you have a contradiction.  The interaction they found was significant.   You can get the four effects they report by adding intercept to each of the RANDOM statements:   proc glimmix data=d method=laplace; class emotion_type target_video_num emotion_type participant_id; model correct_recognized = emotion_type cent_trans_ART emotion_type*cent_trans_ART/dist=binary; random intercept emotion_type/subject=participant_id type=chol; /* and here is the hard part, fitting a continuous variable as a random_effect*/ random intercept cent_trans_ART/subject=target_video_num; /*insert any lsmeans/oddsratio estimate/lsmeans/lsmestimate statements here, or use a STORE statement*/ run; To get the cent_trans_ART variance components, you may have to create an identical variable (say cent_trans_ART_cat =cent_trans_ART) and include it in the class statement and in the second RANDOM statement.  That is going to push the execution time up immensely and may not converge, even in geological time.  There are other ways, but they aren't going to give you identical output.   SteveDenham   ... View more

Hi Eman,  I followed it all, and thought I could write code, but I ran into an issue with emotion_type as a repeated measure.  I know it is repeated within a subject, but is there any ordering to the presentation?  If so, is it identical for all subjects?  Judging from the R model, I would say that there is no ordering, or else there would be an additional R side factor covering the temporal effect.  So, in the R model emotion_type is a G side effect, with a two level unstructured covariance structure. I might be tempted to see if the simpler compound symmetry effect improved the fit based on corrected AIC.  So here is my proposed code:   proc glimmix data=d method=laplace; class emotion_type target_video_num emotion_type participant_id; model correct_recognized = emotion_type cent_trans_ART emotion_type*cent_trans_ART/dist=binary; random emotion_type/subject=participant_id type=chol; /* and here is the hard part, fitting a continuous variable as a random_effect*/ random cent_trans_ART/subject=target_video_num; /*insert any lsmeans/oddsratio estimate/lsmeans/lsmestimate statements here, or use a STORE statement*/ run; This fits an unequal slopes/unequal intercept model for cent_trans_ART.  A reduced model may be appropriate if the emotion_type by cent_trans_ART interaction is not significant.  Oh yeah, no guarantees that this code will converge, etc.   SteveDenham   ... View more

Well, perhaps Hosmer-Lemeshow isn't the best choice for goodness of fit.  Per the documentation:   "Hosmer and Lemeshow (2000) proposed a statistic that they show, through simulation, is distributed as chi-square when there is no replication in any of the subpopulations" I suspect that your categorical variables lead to substantial replication within some of the subpopulations.   When you do have replication and a fair sized sample, check on using the AGGREGATE statement so you can get the Pearson and Deviance goodness of fit tests.  If the data get sparse with all of this, the Details section outlines several other methods.   SteveDenham ... View more

I almost always prefer a modeling approach to questions like this, but this looks to me to be a classic Cochran-Mantel-Haenszel situation.  Stratifying by age category, you get five 2x2 tables (immigrant status by injection status).  Have you looked at this approach?   proc freq data=yourdata; tables agecat*imm_cat*doc_visit/all commonriskdiff; run; I realize this doesn't address gender, income, time since immigration, etc., but this high level analysis will point out whether you need to include other factors.  Since you are employing some sort of stepwise variable selection, it might be good to see first if there is a common risk difference across the age categories.  Then a deeper dive could explain why it differs.   SteveDenham   ... View more