Questions regarding interpretation of linear mixed model

[Junemoe Jeong]

Hello,

I am using Stata/MP 14.1 for Mac. I've been trying to analyze my data from repeated measure design experiment (n=6) with linear mixed models (Command mixed) under the setting like below.

- Dependent variable : volume_after
- Explanatory variable : approach (left, right), pressure (0, 5, 10, 15)
- Group variable : Index

And for better small sample analysis, I used "reml" and "dfmethod (kroger)" option described in Stata material.
(http://www.stata.com/meeting/columbus15/abstracts/materials/columbus15_yang.pdf)

I followed "top-down analysis" described in the textbook (Linear mixed models : A practical guide using statistical software, second edition).

Among three models described below, "lrtest" showed that Model 2. the better fitted model.
Model 1. mixed volume_after i.pressure approach_d i.pressure#c.approach_d || index: , covariance (identity) variance reml dfmethod(kroger)
Model 2. mixed volume_after i.pressure approach_d i.pressure#c.approach_d || index: pressure, covariance (unstruct) variance reml dfmethod(kroger)
Model 3. mixed volume_after i.pressure approach_d i.pressure#c.approach_d || index: pressure, covariance (unstruct) variance reml dfmethod(kroger) residuals(independent, by(approach_d))

* "approach_d" is a dummy variable for "approach"

The result of Model 2

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -69.838378  
Iteration 1:   log restricted-likelihood = -69.684868  
Iteration 2:   log restricted-likelihood = -69.680506  
Iteration 3:   log restricted-likelihood = -69.680503  

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression                   Number of obs     =         20
Group variable: index                           Number of groups  =          6

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.3
                                                              max =          4
DF method: Kenward-Roger                        DF:           min =       4.00
                                                              avg =       4.67
                                                              max =       5.28

                                                F(7,     4.41)    =      25.93
Log restricted-likelihood = -69.680503          Prob > F          =     0.0023

---------------------------------------------------------------------------------------
         volume_after |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
----------------------+----------------------------------------------------------------
             pressure |
                   5  |   90.16037   65.08291     1.39   0.223    -75.89707    256.2178
                  10  |    368.853   121.7059     3.03   0.039     30.94312    706.7628
                  15  |   1015.671   187.8344     5.41   0.005     506.6321    1524.709
                      |
           approach_d |    -3.8654   28.93499    -0.13   0.899    -77.06403    69.33323
                      |
pressure#c.approach_d |
                   5  |   -12.2128   92.04114    -0.13   0.899    -247.0535    222.6279
                  10  |   234.4865   172.1182     1.36   0.245    -243.3902    712.3632
                  15  |   31.96351   265.6379     0.12   0.910    -687.9255    751.8526
                      |
                _cons |   17.02083   20.46012     0.83   0.441    -34.73841    68.78008
---------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
index: Unstructured          |
               var(pressure) |    423.062   313.9038      98.81635    1811.253
                  var(_cons) |   190.4463   433.3242       2.20306    16463.38
         cov(pressure,_cons) |  -283.8496   388.1759      -1044.66    476.9611
-----------------------------+------------------------------------------------
               var(Residual) |   1065.404   532.7942      399.7972    2839.153
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 17.49                 Prob > chi2 = 0.0006

Note: LR test is conservative and provided only for reference.

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |         20         .   -69.6805      12     163.361   175.3098
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

Finally, I want to reduce non-significant fixed effect "approach_d" & "pressure#c.approach_d" with "testparm" command.

. testparm i.pressure#c.approach_d

 ( 1)  [volume_after]5.pressure#c.approach_d = 0
 ( 2)  [volume_after]10.pressure#c.approach_d = 0
 ( 3)  [volume_after]15.pressure#c.approach_d = 0

           chi2(  3) =   28.59
         Prob > chi2 =    0.0000

. testparm approach_d

 ( 1)  [volume_after]approach_d = 0

           chi2(  1) =    0.02
         Prob > chi2 =    0.8937

. testparm i.pressure

 ( 1)  [volume_after]5.pressure = 0
 ( 2)  [volume_after]10.pressure = 0
 ( 3)  [volume_after]15.pressure = 0

           chi2(  3) =  107.60
         Prob > chi2 =    0.0000

Here're questions...

1) With the "testparm" results above, can I just remove "i.pressure#c.approach_d" and get final model for estimates?

Final model(?)
mixed volume_after i.pressure i.pressure#approach_d ||index: pressure, covariance (identity) variance dfmethod(kroger) reml

2) If so, how can I understand my final model? I don't think I can say "approach is not relevant to volume_after".

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -70.546237  
Iteration 1:   log restricted-likelihood = -70.546229  
Iteration 2:   log restricted-likelihood = -70.546229  

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression                   Number of obs     =         20
Group variable: index                           Number of groups  =          6

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.3
                                                              max =          4
DF method: Kenward-Roger                        DF:           min =       4.21
                                                              avg =       6.30
                                                              max =       9.69

                                                F(7,     7.13)    =      29.66
Log restricted-likelihood = -70.546229          Prob > F          =     0.0001

-------------------------------------------------------------------------------------
       volume_after |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
           pressure |
                 5  |   90.16037   58.01232     1.55   0.165     -47.3103     227.631
                10  |    368.853   104.6488     3.52   0.019     92.36201    645.3439
                15  |   1007.622   161.5402     6.24   0.001      594.131    1421.112
                    |
pressure#approach_d |
               0 1  |    -3.8654   32.23463    -0.12   0.907    -76.00165    68.27085
               5 1  |   -16.0782   78.07926    -0.21   0.845    -213.4368    181.2804
              10 1  |   230.6211   145.8365     1.58   0.185    -166.5587    627.8009
              15 1  |   24.77871   227.0594     0.11   0.917    -562.6496     612.207
                    |
              _cons |   17.02083   22.79332     0.75   0.473     -33.9872    68.02887
-------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
index: Identity              |
         var(pressure _cons) |    303.438   195.4811      85.84362    1072.585
-----------------------------+------------------------------------------------
               var(Residual) |   1255.169   657.3175      449.7154    3503.214
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 15.76         Prob >= chibar2 = 0.0000

3) I just tried Model 4. and compared with Final model (?) with AIC, BIC & "lrtest", which showed that final model is better than model 4.
** "lrtest" was performed without the option "reml" & "dfmethod (kroger)"

. lrtest model_final model4

Likelihood-ratio test                                 LR chi2(8)  =     22.32
(Assumption: model4 nested in model_final)                 Prob > chi2 =    0.0044

But, it seems that the result of Model 2. showed "approach_d" and "i.pressure#c.approach_d" is not signifincant"

How can I understand this result?

4)
I thought this problem is due to the small sample size,
I tried "Mann-Whitney U test" to remove "approach" factor before performing linear mixed model....then Model 4 will be my final model.

What are the limitations in this approach compared to the approach above? If I cannot increase my sample size, can I adopt this approach?


Model 4. mixed volume_after i.pressure || index: pressure, covariance (unstruct) variance reml dfmethod(kroger)

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -98.097027  
Iteration 1:   log restricted-likelihood = -97.956087  
Iteration 2:   log restricted-likelihood = -97.946855  
Iteration 3:   log restricted-likelihood = -97.946798  
Iteration 4:   log restricted-likelihood = -97.946784  
Iteration 5:   log restricted-likelihood =  -97.94678  

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression                   Number of obs     =         20
Group variable: index                           Number of groups  =          6

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.3
                                                              max =          4
DF method: Kenward-Roger                        DF:           min =       5.19
                                                              avg =       6.78
                                                              max =       8.83

                                                F(3,     7.27)    =      25.98
Log restricted-likelihood =  -97.94678          Prob > F          =     0.0003

------------------------------------------------------------------------------
volume_after |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    pressure |
          5  |   84.05397   54.81305     1.53   0.160    -40.29687    208.4048
         10  |   486.0962    92.6376     5.25   0.003     250.5715    721.6209
         15  |   1023.594    148.682     6.88   0.000     659.7686     1387.42
             |
       _cons |   15.08813   25.43626     0.59   0.571    -44.90898    75.08525
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
index: Unstructured          |
               var(pressure) |   446.1804   324.3317      107.3411     1854.62
                  var(_cons) |   445.8629   1055.327      4.309887    46125.04
         cov(pressure,_cons) |  -446.0216    638.208     -1696.886    804.8431
-----------------------------+------------------------------------------------
               var(Residual) |   3436.156   1472.698      1483.397    7959.545
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 11.67                 Prob > chi2 = 0.0086

Note: LR test is conservative and provided only for reference.

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |         20         .  -97.94678       8    211.8936   219.8594
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

Thanks,

Jay

[Clyde Schechter]

Before getting into the questions you posed, there are aspects of your models that are peculiar at least, and probably improper specifications. On the one hand, in the fixed-effects you treat pressure as a discrete variable (fair enough if you expect highly non-linear effects of pressure over the four values available). But then you also include random slopes on pressure in models 2, 3, and 4. That's incoherent. You need to treat pressure the same way in both parts of the model. Since the random effects part of the model specification for -mixed- does not support factor variable notation, you need to generate your own indicator variables for the levels of pressure. Actually, probably the simplest way to go about this is to drag out the otherwise largely obsolete -xi:- prefix. So something like -xi: mixed volume_after i.pressure*approach_d || index: i.pressure-. Whether you go this particular route or make your model internally consistent in some other way (such as adding a c.pressure term to the fixed effects while also keeping the i.pressure effects and interactions, or something else entirely) you need to fix this problem to get results that can be interpreted.

It's also strange that you are treating approach_d as a continuous variable by specifying the c. prefix when you interact it with pressure, when it is a dichotomy. This won't lead to any incorrect results in a linear model, but it is strange.

Another general concern is doing this modeling with only 6 index: groups and a total of 20 observations. With 6 index groups you really don't get much sampling at the group level, so the estimates of the random intercepts and slopes will not be very reliable. Moreover, remember that ML estimation is asymptotically valid -- this sample strikes me as too small to really qualify.

Putting aside the issue noted above, below is some advice regarding your questions. My comments are meant as general advice about interpreting these kinds of models. I would not attempt to actually draw conclusions from the models shown in your post, given the limitations I have noted above.

1) With the "testparm" results above, can I just remove "i.pressure#c.approach_d" and get final model for estimates?

Final model(?)
mixed volume_after i.pressure i.pressure#approach_d ||index: pressure, covariance (identity) variance dfmethod(kroger) reml

2) If so, how can I understand my final model? I don't think I can say "approach is not relevant to volume_after".

Probably not. The model including the pressure # approach interaction without including the approach main effect is a model in which you are constraining the effect of approach to be 0 conditional on pressure = 0. If that's really what you want, OK, but it's an odd model, and in any event I don't think it's what you have in mind. Be that as it may, you definitely cannot conclude that approach is not relevant to volume_after because you have significant interaction terms involving approach--look at your -testparm- results shown just after. Also, to really conclude that approach is irrelevant to volume_after you would need to look at -testparm approach_d i.pressure#c.approach_d-, a test that does not appear in your output. If that test were negative, you could consider removing not just approach_d but approach_d and the interaction terms from the model, and that would make some sense. (And it appears that is what you did in Model 4.)

3) I just tried Model 4. and compared with Final model (?) with AIC, BIC & "lrtest", which showed that final model is better than model 4.
** "lrtest" was performed without the option "reml" & "dfmethod (kroger)"
... [code block omitted]
But, it seems that the result of Model 2. showed "approach_d" and "i.pressure#c.approach_d" is not signifincant"

How can I understand this result?

But in fact you did not show that approach_d and i.pressure#c.approach_d are not significant in model 2; the -testparm- results for those show that they are highly jointly significant (even though individually they are not). The LR test results are consistent with that.

[Junemoe Jeong]

Thank you for kind comment, Clyde.

Aside from techinical aspects that I should consider, I need your advise what I can get from this sample with linear mixed model.

I understand that this small size model has statical limitation but, I cannot make this sample size meaningfully increased due to the characteristics of my areas. (most of experiments in my field performed within 10 samples) Also, this research has an exploratory purpose, which essentially has unexpected uncontrolled factors.

Within this limitations, therefore, my initial plan was to remove "approach" effect using "Mann-Whitney U test" before construct a linear mixed model.
And wanted to make a minimal model to evaluate the effects of pressure on the total volume controlling variation of individual subject.

My interest is to find most effective pressure interval.

So, if I want to get contrast to compare pressure intervals, will it be valid to run the model without random slope like below?

. mixed volume_after i.pressure || index: , covariance (identity) variance reml dfmethod(kroger)

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log restricted-likelihood =  -103.6278  
Iteration 1:   log restricted-likelihood = -103.62459  
Iteration 2:   log restricted-likelihood = -103.62459  

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression                   Number of obs     =         20
Group variable: index                           Number of groups  =          6

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.3
                                                              max =          4
DF method: Kenward-Roger                        DF:           min =      11.03
                                                              avg =      12.62
                                                              max =      15.20

                                                F(3,    11.82)    =      32.40
Log restricted-likelihood = -103.62459          Prob > F          =     0.0000

------------------------------------------------------------------------------
volume_after |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    pressure |
          5  |   84.05397   69.94044     1.20   0.255    -69.83682    237.9448
         10  |   486.0962   69.94044     6.95   0.000     332.2054     639.987
         15  |   834.2477    107.181     7.78   0.000     603.0616    1065.434
             |
       _cons |   15.08813    53.7287     0.28   0.783    -99.29776     129.474
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
index: Identity              |
                  var(_cons) |   2645.645   5295.144      52.34654    133713.4
-----------------------------+------------------------------------------------
               var(Residual) |      14675   6382.706      6256.891    34418.94
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 0.32          Prob >= chibar2 = 0.2868

. margins pressure, pwcompare (effects) mcompare (bonferroni)

Pairwise comparisons of adjusted predictions

Expression   : Linear prediction, fixed portion, predict()

---------------------------
             |    Number of
             |  Comparisons
-------------+-------------
    pressure |            6
---------------------------

------------------------------------------------------------------------------
             |            Delta-method    Bonferroni           Bonferroni
             |   Contrast   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    pressure |
    5 vs  0  |   84.05397   69.94044     1.20   1.000    -100.4669    268.5749
   10 vs  0  |   486.0962   69.94044     6.95   0.000     301.5753    670.6171
   15 vs  0  |   834.2477   101.7628     8.20   0.000     565.7713    1102.724
   10 vs  5  |   402.0422   69.94044     5.75   0.000     217.5214    586.5631
   15 vs  5  |   750.1938   101.7628     7.37   0.000     481.7173     1018.67
   15 vs 10  |   348.1515   101.7628     3.42   0.004     79.67509     616.628
------------------------------------------------------------------------------

[Clyde Schechter]

Yes, this model makes sense, but whether you can eliminate approach_d remains in question.

I think a Mann-Whitney isn't a valid method here because the observations are not independent--one of the few assumptions required for this kind of non-parametric test.

Given that random-effects estimation with 6 groups is iffy, I would probably start out with a purely fixed-effects model like:

regress volume_after i.pressure##i.approach_d i.index
testparm 1.approach_d 1.approach_d#i.pressure

If the results of -testparm- are not statistically significant, and if none of the coefficients of those terms looks too large to ignore, then I would take out approach and the interactions with pressure, and go to the model you have in #3, or perhaps to a version treating index as a fixed effect only:

regress volume_after i.pressure i.index

If one of your concerns is that the effects of the different levels of pressure may vary among the 6 replications, it would be tempting to do the fixed-effects version of a random slopes model, namely:

regress volume_after i.pressure##i.index

but if I've done my arithmetic correctly your data set won't give you that many degrees of freedom, so you would have to settle for a true random slopes model such as:

xi: mixed volume_after i.pressure || index: i.pressure // FEEL FREE TO ADD OPTIONS HERE

While this is unsatisfactory due to the small number of group and small sample, I think it's probably the best one can do under the circumstances.

[Junemoe Jeong]

Thank you again, I really appreciate it

Based on your comment, I can take out approach and the interactions with pressure.

Then tried other options you suggested. But, I don't know what criteria can be used for finding the most fitted model.

Model 1. (from Posting #3)
Model 2. only fixed effects
Model 3. + fixed-effects version of a random slopes model
Model 4.1 Without options ; ML estimation (Default)
Model 4.2. REML estimation
Model 4.3.ML estimation + covariance structure (identity)
Model 4.4. ML estimation + covariance structure (unstruct)

What I can think of is that I can choose the one showed p-value < 0.5 in every part of the model.

Results can be summarized :
- Fixed effect & random effect ; p <0.5 :: Model 2. Model 4.1.
- Random effect estimation ; p > 0.5 :: Model 1. ,Model 4.3.
- Cannot get results :: Model 3, Model 4.2. 4.4

(I attached the output at the bottom of this posting)

Therefore, I can find Model 2 & Model 4.1. will be the most fitted model.

And I want to include pressure as a random effect, which makes Model 4.1. is my final model.

But the problem is that I can't get contrast with "margin" command for comparison between pressure intervals.

. margins pressure, pwcompare (effects) mcompare (bonferroni)
factor 'pressure' not found in list of covariates
r(322);

--------------------------------------------- Models Outputs ---------------------------------------------

Model 2.

regress volume_after i.pressure i.index

      Source |       SS           df       MS      Number of obs   =        20
-------------+----------------------------------   F(8, 11)        =     14.03
       Model |  1618227.78         8  202278.472   Prob > F        =    0.0001
    Residual |  158628.238        11  14420.7489   R-squared       =    0.9107
-------------+----------------------------------   Adj R-squared   =    0.8458
       Total |  1776856.01        19  93518.7375   Root MSE        =    120.09

------------------------------------------------------------------------------
volume_after |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    pressure |
          5  |   84.05397   69.33193     1.21   0.251    -68.54458    236.6525
         10  |   486.0962   69.33193     7.01   0.000     333.4977    638.6948
         15  |   879.8903   105.9063     8.31   0.000     646.7921    1112.988
             |
       index |
          2  |   -14.2096   84.91392    -0.17   0.870    -201.1039    172.6847
          3  |   125.8611   94.93665     1.33   0.212    -83.09305    334.8153
          4  |   149.2197   94.93665     1.57   0.144    -59.73448    358.1738
          5  |   153.4674   94.93665     1.62   0.134    -55.48675    362.4216
          7  |  -34.99855   94.93665    -0.37   0.719    -243.9527    173.9556
             |
       _cons |  -48.13521   76.21258    -0.63   0.541     -215.878    119.6075
------------------------------------------------------------------------------

Model 3.

regress volume_after i.pressure##i.index => No results


Model 4.

xi: mixed volume_after i.pressure || index: i.pressure // FEEL FREE TO ADD OPTIONS HERE
Model 4.1) Without options.

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =         20
Group variable: index                           Number of groups  =          6

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.3
                                                              max =          4

                                                Wald chi2(3)      =     865.59
Log likelihood = -98.649616                     Prob > chi2       =     0.0000

-------------------------------------------------------------------------------
 volume_after |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
 _Ipressure_5 |   84.05397   22.84347     3.68   0.000     39.28158    128.8264
_Ipressure_10 |   486.0962   84.46272     5.76   0.000     320.5523    651.6401
_Ipressure_15 |   810.6421   28.32737    28.62   0.000     755.1215    866.1627
        _cons |   15.08813   1.257416    12.00   0.000     12.62364    17.55262
-------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
index: Independent           |
               var(_Ipre~_5) |   3130.946   1807.406      1009.952     9706.22
               var(_Ipre~10) |   42803.71   24712.74      13805.11    132715.9
               var(_Ipre~15) |    1604.88   1604.882      226.0685    11393.18
                  var(_cons) |   9.486566   5.476965      3.059688    29.41311
-----------------------------+------------------------------------------------
               var(Residual) |   1.40e-11   7.20e-11      5.74e-16    3.40e-07
------------------------------------------------------------------------------
LR test vs. linear model: chi2(4) = 50.11                 Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

Model 4.2) With option "reml"
Model 4.4) With covariance (unstruct)
-> Both Error.

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -84.152197  
Iteration 1:   log restricted-likelihood = -83.624036  
Iteration 2:   log restricted-likelihood = -83.610082  
Iteration 3:   log restricted-likelihood = -83.584251  (not concave)
Iteration 4:   log restricted-likelihood = -83.583362  
Iteration 5:   log restricted-likelihood = -83.575105  
Iteration 6:   log restricted-likelihood = -83.552295  
Iteration 7:   log restricted-likelihood = -83.550809  
Iteration 8:   log restricted-likelihood = -83.550154  
Iteration 9:   log restricted-likelihood = -83.550062  
Iteration 10:  log restricted-likelihood = -83.549989  (not concave)
Iteration 11:  log restricted-likelihood = -83.549986  (not concave)
Iteration 12:  log restricted-likelihood = -83.549986  (not concave)
Iteration 13:  log restricted-likelihood = -83.549986  (not concave)
Iteration 14:  log restricted-likelihood = -83.549986  (not concave)
Iteration 15:  log restricted-likelihood = -83.549986  (not concave)
Iteration 16:  log restricted-likelihood = -83.549986  (not concave)
Iteration 17:  log restricted-likelihood = -83.549986  (not concave)
Iteration 18:  log restricted-likelihood = -83.549986  (not concave)
--Break--
r(1);

Model 4.3) With option "covariance (identity)"

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log likelihood = -122.93987  (not concave)
Iteration 1:   log likelihood = -122.78353  
Iteration 2:   log likelihood = -122.70596  
Iteration 3:   log likelihood =  -122.7041  
Iteration 4:   log likelihood = -122.70409  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =         20
Group variable: index                           Number of groups  =          6

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.3
                                                              max =          4

                                                Wald chi2(3)      =     163.97
Log likelihood = -122.70409                     Prob > chi2       =     0.0000

-------------------------------------------------------------------------------
 volume_after |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
 _Ipressure_5 |   84.05397   52.43484     1.60   0.109    -18.71642    186.8244
_Ipressure_10 |   486.0962   52.43484     9.27   0.000     383.3258    588.8666
_Ipressure_15 |    844.875   83.57001    10.11   0.000     681.0808    1008.669
        _cons |   15.08813   42.60227     0.35   0.723    -68.41078    98.58704
-------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
index: Identity              |
     var(_Ipre~_5.._cons)(1) |   5282.966     4776.6      897.9724    31080.84
-----------------------------+------------------------------------------------
               var(Residual) |   5606.753   4867.978      1022.513    30743.56
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 2.00          Prob >= chibar2 = 0.0786
(1) _Ipressure_5 _Ipressure_10 _Ipressure_15 _cons

[Clyde Schechter]

Then tried other options you suggested. But, I don't know what criteria can be used for finding the most fitted model.

Model 1. (from Posting #3)
Model 2. only fixed effects
Model 3. + fixed-effects version of a random slopes model
Model 4.1 Without options ; ML estimation (Default)
Model 4.2. REML estimation
Model 4.3.ML estimation + covariance structure (identity)
Model 4.4. ML estimation + covariance structure (unstruct)

What I can think of is that I can choose the one showed p-value < 0.5 in every part of the model.

I really would not do that. Given 20 observations in 6 groups and 3 parameters being fit in the smallest model, you have a noisy data generating process. The best fitting model is quite likely to simply be an overfitting of the random errors, rather than a model of the process that would replicate were the study to be repeated. Moreover, picking the model with the most significant p-values completely defeats the very meaning of p-values! It's called p-hacking. Although it's widely practiced, most will agree that it's bad practice at best. In some circles it has come to be viewed as scientific misconduct!

Given the very limited size of your data set, there is really no defensible way to select the "best" model. Moreover, what really matters are your estimates of the effects of pressure at each level: and those are the same in all of these models. The rest of the modeling tries to pin down variation in those effects by group--but at the end of the day you just don't have enough data to do that in a serious way. What you can say is that in model 2 in #5 suggest that there isn't that much group (index)-level variation. Without even looking at p-values, you can see in model 2 that the group level effects are small relative to the differences in the pressure level effects, so probably not terribly important. In the model in #3, which is a random-effects version of that, again, the standard deviation of the random effects is sqrt(2645.645) which is approximately 50, and is, again, small relative to the differences in effects of pressure levels. So both approaches suggest that group-level variation isn't playing much of a role here. In a study with more groups one might be able to say that with more force and precision, but I think this is the best we can do with this data set, however strongly we might desire to say something more convincing. So I think the results to focus on are the effects of the different power levels and their magnitudes. The p-values and confidence intervals are simply not reliable in this data structure.

If I were forced at gun point to choose one model, I would go with Model 2 (from #5): not because of the results it gives, but because its a model whose assumptions are least violated by the structure of the data.

Let me comment on what happened with some of the individual models:

Model 3 gave no results because the data set is too small for that many degrees of freedom. I should have done that calculation before recommending you try it: it was perfectly predictable that with 20 observations this model could not be estimated.

Model 4.1: you can't use -margins- after that -xi: mixed- command because pressure is not entered as a factor variable in that model. (When you use the xi: prefix, even though i.pressure looks like a factor variable, it isn't: it's just a shorthand for the varlist _Ipressure*.

Models 4.4 did not converge because the model is not identified. The use of unstructured covariance just added a few extra parameters to be estimated, and that was, again too much for this small data set.

I'm not sure why model 4.2 didn't converge. I can't think of a reason why -reml- estimation would fail here.