You are here

Robust ML estimation in OpenMX

8 posts / 0 new
Last post
jam86hk's picture
Offline
Joined: 05/21/2010 - 03:22
Robust ML estimation in OpenMX

Dear all,

I am a new learner of OpenMX and R. I am interested to know if there is alternative estimation (robust ML such as Satorra-Bentler test statisitc, ADF/WLS by Browne) besides normal-theory ML?

Thank you!

Heining Cham

Mike Cheung's picture
Offline
Joined: 10/08/2009 - 22:37
Dear Heining, ADF can be

Dear Heining,

ADF can be implemented quite easily. The results look reasonable. The other approaches may be more challenging.

Cheers,
Mike

library(MASS)

http://r-forge.r-project.org/projects/sem-additions/

It is used to estimate the fourth order moment matrix.

library(sem.additions)
library(OpenMx)

Generate artificial data based on a one-factor CFA

set.seed(100)
n <- 1000

Population values for 1 factor model

L1 <- matrix(c(1,1.2,1.3,0,0,0,0,0,0,1,1.2,1.3),ncol=2)
P1 <- matrix(c(1,0.3,0.3,1), ncol=2)
E1 <- diag(rep(0.5,6))
Sigma <- L1 %&% P1 + E1 # one factor model
my.df <- data.frame(mvrnorm(n=n, mu=rep(0,6), Sigma=Sigma))

Sample covariance matrix

cov.sample <- cov(my.df)

Sample asymptotic covariance matrix

acov <- adf.wmat(my.df)

Sample covariance matrix

D <- mxMatrix("Full", ncol=6, nrow=6, values=c(cov.sample), name="D")

Factor loadings

L <- mxMatrix("Full", ncol=2, nrow=6, values=c(1,1,1,0,0,0,0,0,0,1,1,1),
free=c(rep(T,3),rep(F,6),rep(T,3)), name="L")

Factor correlation matrix

P <- mxMatrix("Stand", ncol=2, nrow=2, value=0.5,free=TRUE, name="P")

Error covariance matrx

E <- mxMatrix("Diag", ncol=6, nrow=6, values=0.2, free=TRUE, name="E")

Discrepancy between sample covariance matrix and model implied covariance matrix

S <- mxAlgebra( D-(L %&% P + E ), name="S")

Convert it into a vector of covariance elements

vech_S <- mxAlgebra(vech(S), "vech_S")

Weight matrix multiplied by n

acov_S <- mxMatrix("Full", ncol=21, nrow=21, values=n*c(solve(acov)), free=FALSE, name="acov_S")

Objective function

obj <- mxAlgebra( t(vech_S) %&% acov_S , name="obj")
objective <- mxAlgebraObjective('obj')

model <- mxModel("WLS", D, L, P, E, S, vech_S, acov_S, obj, objective, mxCI(c("P")))
fit <- mxRun(model, intervals=TRUE)
summary(fit)

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
Very cool. One small comment

Very cool. One small comment is that in OpenMx 0.3.1 confidence intervals for all the cells of the "P" matrix are calculated, including those cells with fixed values. In OpenMx 0.3.2 (coming soon) only cells with free parameters will be calculated. In the meantime, use mxCI(c("P[1,2]", "P[2,1]")).

I used the following steps to install the sem.additions library. I needed the prerequisites sem, mvnormtest, sna, matrixcalc, and igraph on my system. Then I downloaded the source for sem.additions, opened the tarball, and ran "R CMD INSTALL ." They have binaries for Windows and OS X platforms.

Steve's picture
Offline
Joined: 07/30/2009 - 14:03
Hi Mike, Very cool! This is

Hi Mike,

Very cool! This is the way we hoped OpenMx would be used. Easy implementation of methods as mxAlgebraObjectives and then later incorporating these into optimized C versions in the computational engine.

I hope this is just the first of many discussions that implement well-known methods.

We also hope that people will use the mxAlgebraObjective to experiment with novel methods. If we have accelerated the pace of your development of new methods, the OpenMx team will feel a real sense of accomplishment in our work.

Mike Cheung's picture
Offline
Joined: 10/08/2009 - 22:37
Hi Michael Spiegel and

Hi Michael Spiegel and Steve,

Indeed, I have to thank the core team in developing OpenMx and entertaining requests from end users like me. I asked a few questions related to the vech and vechs calculations in OpenMx. Then vech() and vechs() were implemented a few days later!

Mike

neale's picture
Offline
Joined: 07/31/2009 - 15:14
Very nice going! How about

Very nice going! How about the same approach for binary/ordinal data (preferably allowing for missing values)?

jam86hk's picture
Offline
Joined: 05/21/2010 - 03:22
Thank you! I can install

Thank you! I can install sem.additions. I am trying to start from the WLS syntax in Mike's post.

Mike Cheung's picture
Offline
Joined: 10/08/2009 - 22:37
I am writing a package to

I am writing a package to conduct meta-analysis based on OpenMx. One of the functions is WLS that fits correlation and covariance structure analysis. It also calculates some goodness-of-fit indices. It is available for testing at http://courses.nus.edu.sg/course/psycwlm/internet/metaSEM/index.html When the Jacobian matrix is accessible in a later release of OpenMx, I will implement the robust standard error and robust test statistics.

Thanks for the hard work on OpenMx. It's very flexible and extensible.

Regards,
Mike

library(metaSEM)

Sample correlation matrix

R1 <- matrix(c(1.00, 0.22, 0.24, 0.18,
0.22, 1.00, 0.30, 0.22,
0.24, 0.30, 1.00, 0.24,
0.18, 0.22, 0.24, 1.00), ncol=4, nrow=4)
n <- 1000

Calculate the asymptotic covariance matrix of the sample correlation matrix

acovR <- asyCov(R1, n)

P1: Factor variance

P1 <- mxMatrix("Full", ncol=1, nrow=1, value=1, free=FALSE, name="P1")

L1: Factor loadings

L1 <- mxMatrix("Full", ncol=1, nrow=4, value=c(0.3, 0.4, 0.5, 0.4), free=TRUE, name="L1")

Model implied correlation matrix

Please note that error variances are not involved in correlation structure analysis

impliedR1 <- mxAlgebra(L1 %&% P1, name="impliedR1")

wls.fit1 <- wls(S=R1, acovS=acovR, n=n, impliedS=impliedR1, matrices=c(P1, L1), cor.analysis=TRUE, intervals=TRUE)
summary(wls.fit1)

Call:

wls(S = tssem1.obj$pooledS, acovS = tssem1.obj$acovS, n = tssem1.obj$total.n, impliedS = impliedS, matrices = matrices, cor.analysis = cor.analysis,

Coefficients:

Estimate Std.Error lbound ubound z value Pr(>|z|)

L1[1,1] 0.517058 0.023568 0.470865 0.563252 21.9386 < 2.2e-16 ***

L1[2,1] 0.575070 0.023120 0.529756 0.620385 24.8732 < 2.2e-16 ***

L1[3,1] 0.593858 0.025586 0.543710 0.644005 23.2103 < 2.2e-16 ***

L1[4,2] 0.705137 0.014698 0.676329 0.733944 47.9750 < 2.2e-16 ***

L1[5,2] 0.579020 0.016902 0.545893 0.612147 34.2578 < 2.2e-16 ***

L1[6,2] 0.746402 0.013060 0.720805 0.771999 57.1514 < 2.2e-16 ***

L1[7,2] 0.692214 0.013622 0.665515 0.718913 50.8156 < 2.2e-16 ***

L1[8,3] 0.621997 0.052462 0.519174 0.724820 11.8562 < 2.2e-16 ***

L1[9,3] 0.332554 0.032670 0.268522 0.396586 10.1793 < 2.2e-16 ***

P1[1,2] 0.544762 0.025617 0.494554 0.594969 21.2659 < 2.2e-16 ***

P1[1,3] 0.488475 0.055862 0.378987 0.597964 8.7442 < 2.2e-16 ***

P1[2,3] 0.392610 0.040345 0.313536 0.471684 9.7314 < 2.2e-16 ***

---

Signif. codes: 0 '' 0.001 '' 0.01 '' 0.05 '.' 0.1 ' ' 1

Goodness-of-fit indices:

Value

Sample size 2902.0000

Chi-square of target model 384.9062

DF of target model 24.0000

p value of target model 9.344e-67

Chi-square of independent model 4636.5664

DF of independent model 36.0000

RMSEA 0.0720

SRMR 0.0729

TLI 0.8823

CFI 0.9216

AIC 336.9062

BIC 193.5504