Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.

# Re: st: bifactor rotation

 From Phil Schumm To Statalist Subject Re: st: bifactor rotation Date Tue, 19 Mar 2013 22:19:09 -0500

```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:

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/
```