338793_en-US

liliconcom

MIxed Models/REML results: Do we have documentation how the standard error will be calculated

Using SLS with REML, aka mixed effects, a customer get different results in different packages (R/Python/JMP) for the standard error for the confidence interval. (R/Python/JMP: 0.05/0.08/0.002) this is a factor of 40, which is quite high. Do we have any documentation on how it is exactly calculated (formula). The BLUPS itself as well the estimates for both FE and ME are the sema in the different packages, but the prediction interval and the confidence interval.

 

I thought it might have something to do with a different underlying repeated covariance structure as default or the Kacker-Harville-Correction? 

 

 

Adding an example: https://www.rdocumentation.org/packages/doBy/versions/4.6-3/topics/data-dietox

 

- building model  Weight = Time + Feed + Pig (Random effect)

@jiancao @laural_jmp 

 

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

Ah.  Now I see the JMP Data table.  It didn't come through on the email.  Thanks.

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

The standard error of 0.002 is user-calculated from the saved Conditional Mean CI columns. 

 

 

 

 

 

 

 

Re: Mixed Models/REML results: Do we have documentation how the standard error will be calculated

I ask for the R code. 

Re: Mixed Models/REML results: Do we have documentation how the standard error will be calculated

Hi Martin,

 

I believe the method used to compute the covariance matrix of the variance components is on page 1306 of the following paper:

 

Wolfinger, R. D., Tobias, R. D., and Sall, J. (1994). “Computing Gaussian Likelihoods and Their Derivatives for General Linear Mixed Models.” SIAM Journal on Scientific Computing 15:1294–1310.

 

The SAS documentation for Proc Mixed gives a little more details if that is helpful?

SAS Help Center: Default Output

 

Elizabeth might know of another good reference?

 

Laura

@eclaassen 

 

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

@jiancao ,

I'm not seeing this formula column? I'd be very interested in seeing it, because I don't know how you'd compute it without the degrees of freedom, which aren't straightforward.

According to Martin they were comparing the standard erro...

According to Martin they were comparing the standard errors for the confidence interval. So, after saving the conditional mean CI columns I simply take the difference and then use the Tabulate to get the stand error. Since there are different types of CI I was trying to see which CI they were comparing.

 

I also ran the SAS Proc Mixed but the results are slightly different from JMP's. For example, the variance components estimates are not exactly the same.  

Re: According to Martin they were comparing the standard erro...

Ugh. I'd just written a long reply, and when I hit post, I got an authentication error and the reply was lost. I'll try to recreate...and hope we don't get a duplicate post.

That quantity isn't the SE? The difference = obs + crit*se - (obs - crit*se) = 2*crit*se. Taking the variance of that isn't going to be the standard error we're interested in. It'd be the variance of the width of the confidence intervals, which as far as I know isn't anything 'interesting'. If you get the standard error of the predicted values or the confidence bounds, you get an SE ~0.834, which is more like an average SE because each observation has its own SE. 

The discrepancy between SAS and JMP is minor. The (-2 log)likelihoods are the same, so we're running into either a bit of local optima or some numeric fuzz due to the values in the datatable not being nicely rounded. Even the predicted values vary after the 4th decimal place despite the fixed effects being the same.

/* outpm gives the marginal prediction information (values, se, CI) that corresponds to JMPs Mean Confidence Interval */
proc mixed data=sasuser.dietox nobound cl covtest;
class pig;
model weight= time feed / ddfm=kr2 s outpm=krpredmean;
random pig ;
run;

/* outp gives the conditional prediction information that corresponds to JMPs Conditional Mean CI */
proc mixed data=sasuser.dietox nobound;
class pig;
model weight= time feed / ddfm=satterth s outp=spred;
random pig;
run;

Really, any discrepancy in the predicted value SEs between JMP and Python/R is stemming from the difference in the SE for the group/pig variance as noted in Martin's screenshots. The variance of the variance comes into play. I don't know how Python is calculating it, and it seems odd that all of the other SEs agree and that one doesn't. Again, SAS and JMP essentially agree and JMPs SLS-REML and Mixed agree while using slightly different routines.

Re: Mixed Models/REML results: Do we have documentation how the standard error will be calculated

Thanks, @martindemel . I'll take a look.

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

It looks like any discrepancy in the confidence intervals is coming from the computed standard errors of the variance components.
The Python output has the group/pig variance as 31.836 with SE of 2.125. (From Martin's screenshot)
The R output is:

Random effects:
 Groups   Name        Variance Std.Dev.
 Pig      (Intercept) 31.836   5.642   
 Residual              7.364   2.714   
Number of obs: 789, groups:  Pig, 72

JMP has:

SAS has: 

Covariance Parameter Estimates
Cov Parm	Estimate	Standard Error	Z Value	Pr Z	Alpha	Lower	Upper
Pig	31.8276	5.4936	5.79	<.0001	0.05	21.0604	42.5949
Residual	7.3645	0.3898	18.89	<.0001	0.05	6.6565	8.1924

 

All packages obtain the same loglikelihood.
The parameter estimates are essentially the same across all packages.
The SE for the pig/group variance component is comparable for SAS, JMP and R, but very different in Python.
The SE for the residual variance is comparable for SAS and JMP, but R is very different and I see no SE for that estimate in the Python output screenshot.

Because those SEs of the variance components are needed in the further calculations, we get differences in the prediction confidence intervals. I don't see a way to get the degrees of freedom and/or SEs that R is using in the confidence intervals to directly compare those. We only see the resulting discrepancies in the intervals.

To see the degrees of freedom for various intervals, Custom Tests can be added in JMP. This script includes custom tests equivalent to the Mean Confidence Interval and Conditional Mean CI. You can see the predicted value and SE match what is obtained when saving those columns. The additional information is the DFDen used for the test/CI. That DFDen can be used to find the critical value to use with the SE for the CI.

Fit Model(
 Y( :Weight ),
 Effects( :Feed, :Time ),
 Random Effects( :Pig ),
 NoBounds( 1 ),
 Personality( "Standard Least Squares" ),
 Method( "REML" ),
 Emphasis( "Minimal Report" ),
 Run(
  :Weight << {Summary of Fit( 1 ), Analysis of Variance( 0 ),
  Parameter Estimates( 1 ), Scaled Estimates( 0 ),
  Plot Actual by Predicted( 0 ), Plot Regression( 0 ),
  Plot Residual by Predicted( 0 ), Plot Studentized Residuals( 0 ),
  Plot Effect Leverage( 1 ), Plot Residual by Normal Quantiles( 0 ),
  Custom Test(
   [1 5.200005 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0],
   Label( "obs2 marginal" )
  ), Custom Test(
   [1 5.200005 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0],
   Label( "obs2 conditional" )
  )}
 )
);

I guess this 'explains' where the problem is coming from, but doesn't explain why. I don't know how R or Python are obtaining the SEs of the variance components, but there's clearly a difference there. JMP and SAS procedures are fully documented (I think the JMP doc points to the SAS doc for all the deep statistical inference information). Relevant Doc for Viya LMIXED procedure (same as regular PROC MIXED) 

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

Probably not. While technically not the same thing, usually, the terms get used interchangeably.

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

There is no option to change estimation method in JMP.

JMP does not have the option to output the Hessian matrix (we don't store it directly because it's large, typically), but in JMP 16 in Fit Mixed there are now options to include the final asymptotic covariance matrices in the report.

I apologize for the slow reply. Despite 'subscribing' to this topic, I'm not notified very quickly of new replies. If you need an answer more quickly on anything, feel free to tag me. I think I get notified of tags much faster!

Re: MIxed Models/REML results: Do we have documentation how the standard error will be calculated

Elisabeth, no worries, you replied fast enough:)

Thanks a lot for your investigations. This helps a lot to answer the questions and to understand what happened. Though we may not be able to fully uncover some differences in the results between R, Python and JMP it is much clearer now. 

All the best,

Martin