Re: st: bifactor rotation

[Phil Schumm <pschumm@uchicago.edu>]

On Mar 18, 2013, at 6:53 PM, Airey, David C wrote:
> In a class we learned about bifactor rotation in factor analysis. Seems like a good feature to add to Stata if it is not there.

On Mar 19, 2013, at 7:39 AM, Airey, David C wrote:
> Here was a reference:
> 
> http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3253271/pdf/nihms341216.pdf


I've actually been using bifactor models myself recently, and I agree that they can be very useful.  As others have pointed out, a given bifactor model can be fit with -sem-.  Doing a bifactor rotation as described in the Jennrich and Bentler paper you cited above is complicated by the fact that Stata does not expose a general command for rotation according to an arbitrary criterion.  However, those commands (both for orthogonal and oblique rotations) are there under the hood of -rotatemat-.  If you are willing to use undocumented commands (and you should think carefully about the implications of doing so), then you can use these to perform the bifactor rotations described in the paper above (and in another companion paper by the same authors on oblique bifactor rotations).

First, you need to write a command to compute the rotation criterion and its derivative.  Here, I have simply taken the MATLAB code provided in the paper and translated it into Stata:


    clear all
    program bi_quartimin
            args todo q Gq equal L
            tempname G
            mata: bi_quart("`q'","`G'","`L'")
             if (`todo' == 0)  exit
            matrix `Gq' = `G'
    end

    mata:
    void bi_quart(string scalar q, string scalar G, string scalar L)
    {
        real matrix                     Lt, Lt2, N
        real scalar                     k, p
    
        p = rows(st_matrix(L))
        k = cols(st_matrix(L))
        Lt = st_matrix(L)[.,(2..k)]
        Lt2 = Lt:^2
        N = J(k-1,k-1,1) - I(k-1)
        st_numscalar(q, sum(Lt2:*(Lt2*N)))
        st_matrix(G, (J(p,1,0), 4*(Lt:*(Lt2*N))))
    }
    end


Note the the program -bi_quartimin- needs to have the same signature as is used internally for the built-in criteria which -rotatemat- supports (you can figure this out easily from inspection of the file rotatemat.ado).

Now, to illustrate, let's use the first example from the paper:


    mat M = ( 1.17,  0.78,  0.18 \ ///
              2.08,  0.78, -0.22 \ ///
              1.17,  0.78,  0.18 \ ///
              2.15, -0.62, -0.08 \ ///
              1.23, -0.62,  0.32 \ ///
              2.15, -0.62, -0.08)


To rotate this matrix (orthogonally) according to the bifactor criteria, we use:


    findfile rotatemat.ado
    run `"`r(fn)'"'
    tempname T
    Init `T' M
    Minimize `"GPFortho M, crit(bi_quartimin) colnames("1 2 3")"' `T' "protect(25)"


Note that the two programs we need to use (-Init- and -Minimize-) are subroutines of -rotatemat-, so to make them accessible we need to first execute rotatemat.ado using -do- (or -run-, as I've used here).  Once we've done that, we can make the necessary calls.  Doing so yields


    Trial   1 : min criterion   .2107749
    Trial   2 : min criterion   .2107749
    Trial   3 : min criterion   .2107749
    Trial   4 : min criterion   .2107749
    Trial   5 : min criterion   .2107749
    Trial   6 : min criterion   .2107749
    Trial   7 : min criterion   .2107749
    Trial   8 : min criterion   .2107749
    Trial   9 : min criterion   .0000564
    <snip>

    . mat li r(Lh), format(%2.1f)

    r(Lh)[6,3]
          c1    c2    c3
    r1   1.0   0.0   1.0
    r2   2.0  -0.0   1.0
    r3   1.0   0.0   1.0
    r4   2.0   1.0  -0.0
    r5   1.0   1.0   0.0
    r6   2.0   1.0  -0.0


which, as you can see, is the same result given in the paper.  Note that in this case, the first 8 solutions were local maxima, which is something you need to watch out for.


-- Phil


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/