1
$\begingroup$

When I fit a mixed model in the GNU R sommer package, I receive a zero variance components for a main effect (SP2) but high values for the interaction of this main effect (SP2:year).

mod2mmer <- mmer(yield ~ year + location + location:year + env.block, 
random = ~ 
  SP1 + SP1:location + SP1:year + SP1:location:year + 
  SP2 + SP2:location + SP2:year + SP2:location:year +  
  SP1:SP2 + SP1:SP2:location + SP1:SP2:year + SP1:SP2:year:location, data = dat2)

This seems to be a similar problem as in another thread. However, my scenario is not comparable to a genetic study where I want to investigate epistasis and I also wouldn't know how to apply the solution there to a more complex model as I am using. I receive the following output, where e.g. the variance component for SP2 is 0.00 and the interaction SP2:year is 9.09. Looking at the averages of my SP2 effects accross locations and years (in excel) there is clearly variation due to SP2. How can I code that in a correct, yet straigtforward way, so that the main effect is not somehow "absorbed" by the interaction?

======================================================================================
  Multivariate Linear Mixed Model fit by REML                      
***********************************  sommer 4.0  *********************************** 
  ======================================================================================
  logLik     AIC     BIC Method Converge
Value -34.37105 84.7421 117.964     NR     TRUE
======================================================================================
  Variance-Covariance components:
                                      VarComp VarCompSE   Zratio Constraint
SP1.yield-yield                        0.8341    1.0864  0.76769   Positive
SP1:location.yield-yield               0.0000    1.0353  0.00000   Positive
SP1:year.yield-yield                   0.0000    1.0911  0.00000   Positive
SP1:location:year.yield-yield          1.6378    1.4485  1.13072   Positive
SP2.yield-yield                        0.0000    2.4557  0.00000   Positive
SP2:location.yield-yield               0.5784    1.3973  0.41395   Positive
SP2:year.yield-yield                   9.0936    3.3565  2.70922   Positive
SP2:location:year.yield-yield          2.7536    1.8332  1.50209   Positive
SP1:SP2.yield-yield                    0.0700    0.8982  0.07793   Positive
SP1:SP2:location.yield-yield           0.0000    1.2401  0.00000   Positive
SP1:SP2:year.yield-yield               0.0000    1.2517  0.00000   Positive
SP1:SP2:year:location.yield-yield     -1.4621    1.8861 -0.77519   Positive
units.yield-yield                     14.9121    1.3974 10.67147   Positive
======================================================================================
  Fixed effects:
  Trait              Effect Estimate Std.Error t.value
1 yield         (Intercept)   22.876    0.9832  23.267
2 yield            year2019    1.724    1.2919   1.335
3 yield         locationUst   -8.444    1.0580  -7.982
4 yield env.blockFis-2018-2    3.019    0.7244   4.167
5 yield env.blockFis-2019-1   -5.264    1.4570  -3.613
6 yield env.blockFis-2019-2   -6.931    1.4570  -4.757
7 yield env.blockUst-2018-1   -2.788    0.7217  -3.863
8 yield env.blockUst-2019-1   -2.470    0.7050  -3.504
======================================================================================
  Groups and observations:
  yield
SP1                           8
SP1:location                 16
SP1:year                     16
SP1:location:year            32
SP2                          30
SP2:location                 60
SP2:year                     60
SP2:location:year           120
SP1:SP2                     240
SP1:SP2:location            480
SP1:SP2:year                480
SP1:SP2:year:location       960
======================================================================================

xxxxx Edit xxxxx

I did a model selection approach with the step() function in R and effectively many of the interaction terms are irrelevant. Ending up with a model that looks like this

mod2mmer <- mmer(yield ~ year + location + location:year + env.block, 
        random = ~ 
          SP1 + SP1:location:year + SP2 + SP2:location:year + SP2:year, data = dat2) 

After simplifying the model, the main effects of SP2 are now not zero but still extremely small (e.g. -0.1, 0.2, 0.15). In the raw data, the mean values of SP2 scatter around the grand mean (17.5) with much higher differences with a standard deviation of 2.7 (e.g. 5.1, -2.8 etc.*). If I do a linear fixed effects model with the same model parameters (all fixed then obviously) I receive emmeans that come much closer to these means. The BLUPs that I estimate with the mixed model still don't make sense to me, should they not be at least in the range of these values?

*values are means of different levels of SP2 over years and locations

xxxxx Edit2 xxxxx

My output of a model without any interactions in the random term

mod2mmer <- mmer(yield ~ year + location + location:year + env.block, 
    random = ~ 
      SP1 + SP2, data = dat2)

looks like this

================================================================
           Multivariate Linear Mixed Model fit by REML           
************************  sommer 4.0  ************************ 
================================================================
         logLik      AIC      BIC Method Converge
Value -77.44084 170.8817 204.1035     NR     TRUE
================================================================
Variance-Covariance components:
                        VarComp VarCompSE Zratio Constraint
SP1.yield-yield           1.055    0.8749  1.206   Positive
SP2.yield-yield           5.057    1.7382  2.909   Positive
units.yield-yield        21.900    1.4989 14.611   Positive
================================================================
Fixed effects:
      Trait              Effect Estimate Std.Error t.value
1 yield         (Intercept)   22.783    0.8286  27.495
2 yield            year2019    1.880    0.8704   2.160
3 yield         locationUst   -8.506    0.8825  -9.639
4 yield env.blockFis-2018-2    3.088    0.8743   3.533
5 yield env.blockFis-2019-1   -5.336    1.2283  -4.344
6 yield env.blockFis-2019-2   -7.002    1.2283  -5.701
7 yield env.blockUst-2018-1   -2.733    0.8739  -3.128
8 yield env.blockUst-2019-1   -2.470    0.8544  -2.891
================================================================
Groups and observations:
               yield
SP1                8
SP2               30
================================================================
$\endgroup$
8
  • $\begingroup$ Usually this happens when the random structure is over-specified / overfitted. Why do you have so many random effects ? Has the model actually converged ? I don't know the mmer function, so can you explain how the random parameter is used ? Are they all random intercepts ? $\endgroup$ Commented Nov 12, 2020 at 11:08
  • $\begingroup$ Hello Robert, "Why do you have so many random effects ?" This is the full model, including all possible interactions, since I defined SP1 and SP2 as random variables, all the interactions with them also have to be treated as random. "Has the model actually converged ?" Yes "can you explain how the random parameter is used ? Are they all random intercepts ?" Yes, everything after the 'random' part is random (in lmer you would have to put everything in (1|...). $\endgroup$
    – Pisum
    Commented Nov 12, 2020 at 11:24
  • 1
    $\begingroup$ Then I would say you have an over-specified model. Why do you want so many random effects. It doesn't make sense to me. I would almost never include anything that is a fixed effect as, or as part of, a grouping variable for random intercepts, so if you want year and location to be fixed effects, don't include them in the random part. "since I defined SP1 and SP2 as random variables, all the interactions with them also have to be treated as random" - this is a mistake, don't put the interactions in the random part, you are just asking for trouble., $\endgroup$ Commented Nov 12, 2020 at 11:39
  • $\begingroup$ Hello Robert, thanks for your reply, see my edit in the question. $\endgroup$
    – Pisum
    Commented Nov 13, 2020 at 9:44
  • $\begingroup$ Please don't use stepwise procedures. They are very bad. Common sense tells us that those interactions are unnecesssary. It is lucky that the stepwise procedure came the the same conclusion. What does it mean "main effects of SP2 are now not zero" ? SP2 is not a fixed effect. It is a random effect so it makes no sense to talk about the main effect of SP2. $\endgroup$ Commented Nov 13, 2020 at 10:00

1 Answer 1

1
$\begingroup$

as mentioned by Robert the first issue you faced was an overspecified model. When you remove so many of the interaction effects (you could have used also likelihood ratio tests to identify which vc should be removed) you can see that the variation gets absorbed by the different terms in what would be their REML estimates. There's nothing wrong with that, just because you expect the main effect to be more relevant doesn't mean that's the case in a model with a different variance structure, you cannot compare that with the raw data. You fit variance structure models (main effects, diagonal, compound symmetry, unstructured, etc.) depending on your hypothesis and you get REML estimates for those structures so keep that in mind (I recommend you get yourself familiar with variance structures to wrap your mind around it).

Given said that, if you decide you want to force the initial variance components from the main effect, you can fit a second model where you force some initial values from the first model and see if the interactions still explain certain variation. Since you didn't provide your data I will make my point with a genetics example included in the package where you fit a second term but because of the known fact that the second term is not orthogonal to the first one and you expect will absorbe some or all the variation you fixed the variance component from a first model so those don't change using the Gti and Gtc arguments.

data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
#### create the variance-covariance matrix
A <- A.mat(GT) # additive relationship matrix
#### look at the data and fit the model
mix1 <- mmer(Yield~1,
             random=~vs(id,Gu=A),
             rcov=~units,
             data=DT, verbose = FALSE)
summary(mix1)$varcomp
#                     VarComp VarCompSE    Zratio Constraint
# u:id.Yield-Yield   650.4145  325.5562  1.997856   Positive
# units.Yield-Yield 4031.0153  344.6051 11.697493   Positive
####=========================================####
#### adding dominance and forcing the other VC's
####=========================================####
DT$idd <- DT$id;
D <- D.mat(GT) # dominance relationship matrix
mm <- matrix(3,1,1) ## matrix to fix the var comp
mix2 <- mmer(Yield~1,
             random=~vs(id, Gu=A, Gti=mix1$sigma_scaled$`u:id`, Gtc=mm)
                     + vs(idd, Gu=D),
             rcov=~vs(units),
data=DT, verbose = FALSE)
summary(mix2)$varcomp
#                       VarComp VarCompSE     Zratio Constraint
# u:id.Yield-Yield     650.4145  507.0538  1.2827326      Fixed
# u:idd.Yield-Yield    156.1547  293.1392  0.5326982   Positive
# u:units.Yield-Yield 3946.8524  360.4550 10.9496404   Positive

In you case would be like fitting first a model with SP1 and SP2 and fitting a 2nd model adding the interactions but fixing the ones for SP1 and SP2. Also, once you obtain adjusted means using the predict() function (recommend you to get the newest version in GitHub 4.1.2) whether the variance goes to SP1 or SP1:location the adjusted means for SP1 will include the BLUPs for both terms so it doesn't matter that much where the variation went. All the best.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.