nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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

martindemel

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 @Laura_Lancaster

dietox.jmp

Accepted answers

eclaassen

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.

eclaassen

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)

eclaassen

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!

All comments

Laura_Lancaster

Hi Martin,

Are you able to share the data and the script for the model? That would be super helpful if you can. Also, which standard errors and confidence intervals are different? I'm also attempting to add Elizabeth to this track.

Thanks.

Laura

@eclaassen

Laura_Lancaster

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

eclaassen

Martin,
If R code producing the confidence intervals and standard errors can be provided, that would be helpful. I'm lagging behind in my R knowledge, and have only successfully been able to build the model (using the dietox data and lmer() ) and get the marginal and conditional predicted values. Those values match with JMP.

Using the dietox data, I ran the model specified in both JMP and SAS (PROC MIXED), and the marginal and conditional predicted values and confidence intervals match (when the corresponding DDFM method is specified in SAS - KR2 for marginal and Satterthwaite for conditional). SAS does not have Individual CIs. There was an issue with the degrees of freedom for the Individual CIs that has been fixed in 16.0.

I know the Kackar-Harville corrected covariance matrix is used for tests/intervals for fixed effects only estimates. The Satterthwaite method is used if a random effect is involved.

jiancao

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

martindemel

I had added the screeshots originally, but it did not appear in the post. Here I try it again:

2.125 as Stderr for the variance components

martindemel

I ask for the R code.

martindemel

Regarding the degrees of freedom, we just had this thought, and if they use JMP 15, could this be the reason for the different result? When I run this example in both 16 and 15 I see exactly the same results.

Laura_Lancaster

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

eclaassen

@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.

jiancao

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.

eclaassen

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

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" )
		)}
	)
);

jiancao

I noticed that JMP and SAS label it as Standard Error, while R package calls Std.Dev. Is there a difference? In its fixed effect report R labels as Std. Error.

eclaassen

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

martindemel

Thanks Elisabeth!

Regarding the optimization algorithm used, in the doc it's mentioned Newton-Raphson, but with an otion to be changed in LMIXED. I guess in JMP we just use NR and do not have a change to switch, correct?

And there will also not be any option to write out the Hessian matrix. Having this output we might be able to compare if there has been a difference from itereation to iteration which differs to the R algorithm, allthough there should not be one mathematically. The comparison work I would let the customer do).

Thanks again for your investigations! That was very helpful.