How can I control for convergence in Proc Glimmix procedure for linear mixed models?

[Lyson]

I am working on a linear mixed model to compare results in SAS proc glimmix (or equivalently proc mixed) and CRAN R lme4 package. My challenge is: How do I manipulate the SAS code to get similar (approximately equal) F-ratios for the analysis of variance in SAS and R?

I tried to control for convergence in SAS (see code), but I need help to finetune that SAS code to get the degrees of freedom and/or F-ratio that match those i get in R. 

Data extract (1st 5 observations):

Here is my SAS code:

proc glimmix data=BarleyGroups outdesign=XZ  ;
class factorA factorB  groupA groupB switchA switchB  ;
model yield=  groupA|groupB / ddfm=kr;
random groupA*factorA*switchA; 
random groupA*groupB*factorB*switchB; 
random groupB*factorB*switchB;
random groupB*groupA*factorA*switchA;
random groupA*factorA*switchA*groupB*factorB*switchB;
parms (0.0000) (0.0000) (0.0000) (0.0000) (0.07074)  (0.41458);
nloptions maxiter=1000 maxfunc=10000  xconv=1e-5 gconv=1e-5 technique=nmsimp;
lsmeans  groupA|groupB / lines ;
Title "Linear Mixed Model with Convergence";
run;

Here is my SAS output from the code above:

**********************************************************************************************************************

Here is my R code: 

# Load necessary libraries
library(lme4)
library(emmeans)
library(lmerTest) # Provides Type III tests for fixed effects

control_settings <- lmerControl(
optimizer = "Nelder_Mead", # Use Nelder-Mead optimizer
optCtrl = list(maxfun = 10000), # Set maximum number of function evaluations
check.conv.grad = .makeCC(action = "warning", tol = 1e-5), # Gradient convergence check
check.conv.singular = .makeCC(action = "warning", tol = 1e-5) # Singular fit check
)
Full.mod<-lmer(Yield~groupA+ groupB+groupA:groupB+
(1 | groupA:factorA:switchA) +
(1 | groupB:factorB:switchB) +
(1 | groupA:groupB:factorB:switchB) +
(1 | groupB:groupA:factorA:switchA)+
(1|groupA:factorA:switchA:groupB:factorB:switchB),
data = BarleyGroups,
control = control_settings)
summary(Full.mod)
anova(Full.mod, type=3, ddf = "Kenward-Roger")

Here is the R output:

> anova(Full.mod, type=3, ddf = "Kenward-Roger")
Type III Analysis of Variance Table with Kenward-Roger's method
               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
groupA        0.08694 0.04347     2 4.5808  0.1040 0.9033
groupB        1.51768 0.50589     3 5.5396  1.2056 0.3905
groupA:groupB 1.38705 0.23118     6 7.9107  0.5506 0.7585

************************************************************************************************************************************** 

[jiltao · Accepted solution]

The fact that the standard errors for the covariance parameters are missing indicates that there are convergence issues with your program. I would take care of that first before trying to compare the results with anything else. 

You might start with simplifying your model first. For example, start with simpler random effects. Also try different values for the PARMS statement if necessary. Your NLOPTIONS statement seems to be relaxing the convergence criteria too much. I would try not to have that statement to begin with.

Good luck,

Jill

[Lyson]

Hi Jiltao

Thanks for the response and advice, I am busy trying that one. Convergence
is what I really want to achieve, with reasonable standard errors
reflecting on the output. My model is trying to combine all random effects
into one model, that is the main objective. So I am not quite sure if there
is another approach of forcing this model to behave well without
necessarily simplifying it (excluding other random effects).