Profile likelihood of central estimates with Mixed

[ptadgerv]

Dear SAS community:

Exploring the Profile Likelihood in different procedures covariance parameters is straight forward in NLMIXED, MIXED, and also GLIMMIX. AS example I put this two codes (one for a Normal-normal model in PROC MIXED and the other for a Binomial-Normal model in NLMIXED) of a Random intercept Mixed model, to explore the values for the covariance parameters:

proc nlmixed data=DATA;
PARMS  covsesp=-2 to 2 by .004;
bounds s2usens>=0;bounds s2uspec>=0;
logitp = (msens + usens)*sens + (mspec + uspec)*spec;
p = exp(logitp)/(1+exp(logitp));
model true ~ binomial(n,p);
random usens uspec ~ normal([0 , 0],[s2usens,covsesp,s2uspec])
subject=study_id;run;

proc mixed data=bi_meta NOPROFILE;
class study_id  ;
model logit = dis  non_dis / noint cl df=1000, 1000, 1000, 1000, 1000, 1000;
random dis non_dis / subject=study_id type=un V VCorr;  
repeated / group=rec;
parms  (0.01 to 25.00 by 0.01) (-0.7498) (0.3549)
( 2.11764705882353 ) ( 0.735632183908046 ) ( 0.222222222222222 ) ( 0.380952380952381 ) ( 0.193939393939394 ) ( 2.13333333333333 ) ( 0.094758064516129 ) ( 0.357723577235772 ) ( 0.154411764705882 ) ( 0.392857142857143 ) / eqcons=2 to &nCov_par;
run;

Even further in NLMIXED the PL of the central values or the fixed elements can be explored like this:

proc NLMIXED data=DATA;
PARMS  msens=-1 to 10 by 0.011;
bounds s2usens>=0;bounds s2uspec>=0;
logitp = (msens + usens)*sens + (mspec + uspec)*spec;
p = exp(logitp)/(1+exp(logitp));
model true ~ binomial(n,p);
random usens uspec ~ normal([0 , 0],[s2usens,covsesp,s2uspec])
subject=study_id;run;

My question is:

Can we do the same in PROC MIXED? or the probabilities explorations outside NLMIXED are restricted to the covariance parameters?

[SteveDenham]

I see several approaches in GLIMMIX, but none leap out at me for MIXED. I suppose you could put the estimates into a PARMSDATA dataset and use that option in the PARMS statement (which is the equivalent of the PROC MIXED code you presented).  It just isn't apparent to me what the output you want will be from MIXED.  I really don't see a way to set fixed effects at certain values and look at the likelihoods within MIXED.  The best I can come up with is to generate a range of response variables as Ynew = Yhave - X'beta, where you specify the beta values you would like to profile over, and then fit multiple runs (one for each new set of Ynew).  Does that make sense?

SteveDenham

[ptadgerv]

Dear Steve, 

Thanks for your time and response. I will give more detail to the problem (from a theoric point of view and also from the code)

In this plot, I have been comparing five parameters in different models, mainly Random intercept Mixed models with distributions Normal-Normal, Binomial-Normal, and Bayesian Binomial-Normal and Binomial_mixture of Normal. The 3 first rows are covariance parameters, the last two rows my response variable a logit of a probability in a multivariate hierarchical model. 
The blue lines are the estimates, the curves are the profile likelihood (PL). So far I have been using the book of Russell B. Millar - Maximum Likelihood Estimation and Inference_ With Examples in R, SAS, and ADMB as a guide in the coding process. Also the manuals for MIXED, GLIMMIX and NLMIXED. 

In Proc Mixed, and Glimmix is possible to do a PL of the covariates parameters with this approach (in orange):

ods output Parameters=all21_covsesp;
proc nlmixed data=ma_all cov;
PARMS covsesp=-2 to 2 by .004;
bounds s2usens>=0;bounds s2uspec>=0;
logitp = (msens + usens)*sens + (mspec + uspec)*spec;
p = exp(logitp)/(1+exp(logitp));
model true ~ binomial(n,p);
random usens uspec ~ normal([0 , 0],[s2usens,covsesp,s2uspec])
subject=study_id;run;

data tdata2;do i = 1 to 501;output;end;run;
data tdata2;set tdata2; if i=1 then do;      covp2=-2;   end;   else covp2 + 0.008; drop i;run;

ods output Parameters=gl.all21_covsesp;
PROC GLIMMIX data=ma_all;
class study_id status;
model true/n = status / noint s cl corrb covb ddfm=bw;
random status / subject=study_id S type=chol G; 
covtest tdata=tdata2 / parms;
ods output  covtests=ct;
run;

proc mixed data=bi_meta NOPROFILE;
class study_id  ;
model logit = dis  non_dis / noint cl df=1000, 1000, 1000, 1000, 1000, 1000;
random dis non_dis / subject=study_id V VCorr;  
repeated / group=rec;
parms   (0.9828) (-0.7498) (0.01 to 25.00 by 0.01);
ods output ParmSearch=pl_variance2;
run;

So far all is easy to follow in manuals and the referred book. The issue now is when I want to explore the probability in the average of my responses variable (logit of the probabilities: msens and mspec  in the case of NLMIXED). NLMIXED allows me to provide starting values for my fixed estimations. Like this:

ods output Parameters=all21_msens;
proc nlmixed data=ma_all cov ecov QPOINTS=20 TECH=TRUREG;
PARMS  msens=-1 to 10 by 0.011;
etc...

So far with MIXED or GLIMMIX I only see ways to provide starting values for the covariance parameter, not for the fixed estimation in my hierarchical model. 

Of course, I can pass the Normal-Normal made in MIXED to NLMIXED, but I was trying to avoid this approach. 

Thanks in advance

[SteveDenham · Accepted solution]

Given the output you presented, I can only say to try the following.  You can fix the fixed effects in both MIXED and GLIMMIX by using the workaround I proposed earlier of subtracting x'beta from the observed value, where you set the coefficients in beta to generate your profiles.  I suspect you will need to fit 20 or more of these datasets to get a decent plot.  The MODEL statement would be for an intercept only model (e.g. MODEL ynew=;, with no fixed effects specified). You probably do not want to profile the random effects at the same time (well, perhaps you do, but that is another step), so you will need to add a noiter option to the random statement, and input the values you wish to do the profiling at.  If you do want both at the same time, you'll get a higher dimensional response surface. Use ODS OUTPUT to save the final log likelihood for each of the runs, reshape the data for plotting, and you should have something close to what you want.

SteveDenham