Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Can MIXED estimate an R-only model with heterogeneous compound symmetry?

    I am trying to reproduce the results in Tables 4.1 and 4.3 in Hoffman's (2015) book, Longitudinal Analysis. Her own Stata code (and output) for Chapter 4 can be viewed here, and the data file can be obtained here. The following comment is prominent in her code:

    Code:
    *******   NOTE: NOT ALL MODELS WILL BE POSSIBLE TO ESTIMATE IN STATA     *******
    I believe that one of the models Hoffman could not estimate using Stata is the Compound Symmetry Heterogeneous (CSH) model.

    I am able to reproduce the R and RCORR matrices for Hoffman's simple CS model using this code:

    Code:
    * 2. R-only, compound symmetry (CS) -- or exchangeable in Stata!
    quietly mixed posmood, ///
       || personid: , noconstant reml covariance(unstructured) ///
          residuals(exchangeable,t(studyday))
    estimates store CS
    estat ic, n(200)
    estat wcorrelation, covariance
    estat wcorrelation
    On another page, Hoffman shows this SPSS code to estimate the CSH model she wants:

    Code:
    ECHO 'Ch 4: Empty Means, Compound Symmetry Heterogeneous R-Only Model'.
    MIXED posmood BY PersonID studyday
         /METHOD   = REML
         /PRINT    = SOLUTION TESTCOV R
         /FIXED    =  
         /REPEATED = studyday | COVTYPE(CSH) SUBJECT(PersonID).
    That example, and the following excerpt from the documentation for MIXED (see p. 8) make me wonder if there is a way to estimate the CSH model using by(varname).

    resopts are by(varname) and t(varname).

    by(varname) is for use within the residuals() option and specifies that a set of distinct
    residual-error parameters be estimated for each level of varname. In other words, you use
    by() to model heteroskedasticity
    .
    I wondered if something like this would give me the CSH model:

    Code:
    * 3. Compound symmetry heterogeneous (CSH)
    mixed posmood, ///
       || personid: , noconstant reml covariance(unstructured) ///
          residuals(exchangeable, t(studyday) by(personid))
    But when I execute that code, I see this initial bit of output:

    Code:
    . mixed posmood, ///
    >    || personid: , noconstant variance reml covariance(unstructured) ///
    >       residuals(exchangeable, t(studyday) by(personid))
    Note: t() not required for this residual structure; ignored
    
    Obtaining starting values by EM ...
    
    Performing gradient-based optimization:
    But then the wheels just keep spinning for so long that it seems to me something must be wrong.

    It is entirely possible that I am not using by(varname) correctly. But thus far, I have been unable to find any worked examples (in the documentation or elsewhere) demonstrating its use. Any pointers would be appreciated, therefore.

    I hope I have explained things clearly enough. Thanks for any help you can offer, even if it is only to confirm that Stata's -mixed- command cannot estimate this model.


    PS- I suspect the same issue will come up when I get to a heterogeneous Toeplitz covariance structure later in Table 4.1! Here is the SPSS code for that model:

    Code:
    ECHO 'Ch 4: Empty Means, n-1 Lag Toeplitz Heterogeneous R-Only Model'.
    MIXED posmood BY PersonID studyday
         /METHOD   = REML
         /PRINT    = SOLUTION TESTCOV R
         /FIXED    =  
         /REPEATED = studyday | COVTYPE(TPH) SUBJECT(PersonID).
    --
    Bruce Weaver
    Email: [email protected]
    Version: Stata/MP 18.5 (Windows)

  • #2
    Originally posted by Bruce Weaver View Post
    I believe that one of the models Hoffman could not estimate using Stata is the Compound Symmetry Heterogeneous (CSH) model.

    . . . Thanks for any help you can offer, even if it is only to confirm that Stata's -mixed- command cannot estimate this model.
    You would use this syntax.
    Code:
    mixed posmood ///
       || personid: , reml dfmethod(satterthwaite) ///
          residuals(independent, by(studyday))
    The random effect intercept for person ID gives you the common covariance in the residuals between study days, and the study days' individual independent residual variances give you the heterogeneous (extra) variance over and above the common covariance.

    In order to get the total variance for each study day you'll need to sum the common covariance (the person ID random effects intercept variance) and each study day's independent (individual) residual variance.

    You might find some slight difference from what SPSS or SAS gives you, due to, for example, algorithmic differences in the iterative maximum likelihood method, difference choices as to which convergence criterion governs or differences in default tolerance settings in the convergence criteria. But they should be close in log-likelihood, fixed effect estimates and residual variances and common covariance.

    By the way, I'm not sure why you included covariance(unstructured) for the random effects in the syntax for your two models. When you specify noconstant for the random effects equation and do not specify any random slopes, the entire random effects equation is empty and the covariance() option is meaningless. I recommend omitting it for clarity.

    Comment


    • #3
      Thanks Joseph. If that is what is required using -mixed-, I think I'll just use (SPSS) MIXED if I ever find myself wanting that covariance structure, because it will give me the desired results directly. Ditto for Hoffman's ARH1 and heterogeneous Toeplitz covariance structures, I imagine.
      --
      Bruce Weaver
      Email: [email protected]
      Version: Stata/MP 18.5 (Windows)

      Comment


      • #4
        Originally posted by Bruce Weaver View Post
        . . . I think I'll just use (SPSS) MIXED if I ever find myself wanting that covariance structure . . .
        I think that that's the crux with a lot of these alternative residual covariance structures. I work a biological field, and I suppose that I could come up with a natural phenomenon in that field whose errors are credibly distributed as Compound Symmetry Heterogeneous (CSH) if I thought hard about it. But offhand, I can't think of any.

        I saw your request in the wish list and up-voted it; I agree that Stata ought to match SAS and SPSS if even only for competitive reasons,* but I have to confess that if I have any reason to doubt the adequacy for the task of a simpler residual covariance pattern, then I tend to go straight to the unstructured, or else end up fitting the model using GEE with independent residuals and (cluster) robust standard errors.

        *There is some potential in exploiting the ability to add constraints to the residual variance coefficients in -sem-, but you don't have REML there.

        Comment


        • #5
          Hi Joseph. Thanks for the up-vote in the Wish List thread. To be clear, I would prefer to use Stata -mixed- rather than (SPSS) MIXED. One important reason is that SPSS has no equivalent of -lrtest- for comparing nested models. I agree with your comment about matching SAS & SPSS for competitive reasons.
          --
          Bruce Weaver
          Email: [email protected]
          Version: Stata/MP 18.5 (Windows)

          Comment


          • #6
            Yes, it would be good to have the heterogeneous covariance structures mentioned above in Stata. These are available in R NLME gls() and also in open pharma MMRM packages.

            Heterogeneous Variance: Covariance Structures for Repeated Measures
            Author(s): Russell D. Wolfinger
            Source: Journal of Agricultural, Biological, and Environmental Statistics, Vol. 1, No. 2 (Jun.,1996), pp. 205-230
            Published by: International Biometric SocietyStable
            URL: http://www.jstor.org/stable/1400366 .

            Comment


            • #7
              The CS(H) structure can be estimated by -mixed-, though perhaps the dataset used had idiosyncrasies that prevented the model from converging.

              Code:
              mixed ... , resid(exchangeable, by(group))
              The "time" option isn't needed (-t()-) within residual() because exchangeable implies this information is redundant.

              Comment


              • #8
                Here's a completely silly example to demonstrate the results. Though overall, I agree with supporting more esoteric residual covariance structures. However, I'd much soon have Stata devote their resources to faster estimation routines for all mixed models.

                Code:
                clear *
                cls
                
                local n = 500
                
                mat Means0 = (2,3,3)
                mat Means1 = (2,4,6)
                mat V0 = (1, 0.5, 0.5 \ ///
                          0.5, 1, 0.5 \ ///
                          0.5, 0.5, 1)
                mat V1 = 2 * (1, 0.5, 0.5 \ ///
                              0.5, 1, 0.5 \ ///
                              0.5, 0.5, 1)
                
                mkf X0
                cwf X0
                drawnorm y0 y1 y2, n(`n') double mean(Means0) cov(V0)
                gen byte trt=0
                
                mkf X1
                cwf X1
                drawnorm y0 y1 y2, n(`n') double mean(Means1) cov(V1)
                gen byte trt=1
                
                mkf Data
                cwf Data
                frameappend X0
                frameappend X1
                sort trt, stable
                gen `c(obs_t)' pid = _n
                
                reshape long y , i(trt pid) j(tm)
                
                mixed y i.tm##i.trt || pid : , nocons reml dfmethod(kr) resid(exch, by(trt)) nolog
                Result

                Code:
                . mixed y i.tm##i.trt || pid : , nocons reml dfmethod(kr) resid(exch, by(trt)) nolog
                
                Mixed-effects REML regression                      Number of obs    =    3,000
                Group variable: pid                                Number of groups =    1,000
                                                                   Obs per group:
                                                                                min =        3
                                                                                avg =      3.0
                                                                                max =        3
                DF method: Kenward–Roger                           DF:          min =   995.26
                                                                                avg = 1,383.27
                                                                                max = 1,796.43
                                                                   F(5, 1797.42)    =  1012.34
                Log restricted-likelihood = -4393.3083             Prob > F         =   0.0000
                
                ------------------------------------------------------------------------------
                           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
                -------------+----------------------------------------------------------------
                          tm |
                          1  |   .9340488    .043696    21.38   0.000     .8483022    1.019795
                          2  |   .9379808    .043696    21.47   0.000     .8522342    1.023727
                             |
                       1.trt |  -.0498535   .0772227    -0.65   0.519    -.2013141    .1016071
                             |
                      tm#trt |
                        1 1  |   .9760273    .075682    12.90   0.000     .8275933    1.124461
                        2 1  |   2.996023    .075682    39.59   0.000     2.847589    3.144457
                             |
                       _cons |   2.071528   .0437863    47.31   0.000     1.985604    2.157452
                ------------------------------------------------------------------------------
                
                ------------------------------------------------------------------------------
                  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
                -----------------------------+------------------------------------------------
                pid:                 (empty) |
                -----------------------------+------------------------------------------------
                Residual: Exchangeable,      |
                    by trt                   |
                                   0: var(e) |   .9586191   .0429727      .8779882    1.046655
                                      cov(e) |   .4812838   .0411637      .4006045    .5619631
                                             |
                                   1: var(e) |   2.023054   .0922943      1.850012    2.212282
                                      cov(e) |   1.068447   .0889354      .8941371    1.242758
                ------------------------------------------------------------------------------
                LR test vs. linear model: chi2(3) = 942.89                Prob > chi2 = 0.0000
                Note that you could include a random intercept for the subject, but I strongly recommend against this. Although you end up with the same model, the intercept variance is aliased with the respective variances of the R matrix. If you add these two variances together, you get the same point estimate of total variance, however the standard errors are very badly biased.

                Comment


                • #9
                  Leonardo Guizzetti, isn't there something different between the models under discussion? In one case the different residual variances are with respect to the between group factor while in the other they are with respect to the within group factor?

                  So for example we might have rats measured in two different treatments (treated, control) with independent replicate rats per treatment, and some continuous result is measured in multiple brain regions per rat (repeated measures), and the contrast of interest is treated vs control within brain regions. Variance differs between brain regions, not so much between treatments say.

                  Comment


                  • #10
                    Ah, yes you're right. I hadn't noticed it was a within-subject difference that was requested. In that case I don't see that Stata is capable of this type of CSH structure.

                    Comment


                    • #11
                      Originally posted by Dave Airey View Post
                      . . . some continuous result is measured in multiple brain regions per rat (repeated measures) . . . Variance differs between brain regions, not so much between treatments say.
                      Originally posted by Leonardo Guizzetti View Post
                      In that case I don't see that Stata is capable of this type of CSH structure.
                      Isn't that just this:
                      Code:
                      version 18.0
                      
                      clear *
                      
                      // seedem
                      set seed 1277806009
                      
                      * Rats
                      quietly set obs 32
                      generate byte rid = _n
                      generate double rid_u = rnormal()
                      
                      * Treatments
                      generate byte trt = mod(_n, 2)
                      
                      * Brain regions
                      quietly expand 3
                      bysort rid: generate byte brg = _n
                      
                      * Outcome
                      generate double out = ///
                          0 + rid_u + ///
                          0 * trt + 0 * brg + (trt - 0.5) * (brg - 2) + ///
                              rnormal(0, sqrt(brg))
                      
                      *
                      * Begin here
                      *
                      mixed out i.trt##i.brg || rid: , residuals(independent, by(brg)) ///
                          reml dfmethod(satterthwaite) nolrtest nolog
                      
                      exit
                      Which is the same as that shown above in #2.

                      Comment


                      • #12
                        Thank you for the example, Joseph. Unless I'm mistaken (again), I don't think this is quite the same structure. Following your example code, the model implied covariance and correlation (R and RCORR) are:

                        Code:
                        . estat wcorrelation
                        
                        [output omitted]
                        
                        Correlations:
                        
                                 obs |      1       2       3
                        -------------+-----------------------
                                   1 |  1.000                 
                                   2 |  0.305   1.000         
                                   3 |  0.260   0.200   1.000
                        
                        . estat wcorrelation, covariance
                        
                        Covariances for rid = 1:
                        
                                 obs |      1       2       3
                        -------------+-----------------------
                                   1 |  1.560                 
                                   2 |  0.619   2.642         
                                   3 |  0.619   0.619   3.636
                        Based on the PROC MIXED documentation:
                        - a CS structure, the each variance is the same, and each covariance the same, implying the same correlation at repeated times
                        - a CSH structure, variances are free, and covariance varies but the implied correlation is still constant between repeated times

                        The model above has constant covariance but not correlation.

                        Comment


                        • #13
                          Originally posted by Leonardo Guizzetti View Post
                          Unless I'm mistaken . . .
                          Nope, it is I who stands corrected here. Thanks for the catch.

                          The CSH structure does seem implausible, though: try imagining a data-generating process in nature that has to allow for different variances yet somehow manages to impose a uniform correlation coefficient.

                          Comment


                          • #14
                            I agree with you, Joseph Coveney that the within-subject heterogeneity of CSH is unlikely to arise naturally. One would have to to measure some process over time in equally spaced intervals whose variance increases linearly with time. Off hand I can't think of such a process without some kind of arbitrary selection mechanism to bias the results according to that pattern.

                            The only sensible times I've seen CS used was to assume heterogeneity between groups (as in #8), but I mostly work with biological and medical data.

                            Comment

                            Working...
                            X