[gov2001-l] Problems with Optim

This looks like a rather unorthodox model and without having seen the math it is difficult to assess what exactly the code is trying to accomplish. unfortunately, i also don't know anything about twin studies. a few points: 1. Try the likelihood function first on toy data that you set up so that you know the solution. Only take the code to your real data if you verified your code on a small toy example first 2. In the toy example, try to simplify the problem first, for example try only to maximize using D not M and D jointly. D and M are not parameters but data so you should pass them as an argument directly optim(par.svals,D=rbind(Yd,Xd),M=rbind(Ym,Xm)) 3. triple check your code, for example: vcvd <- matrix(c(sa+sc+se,sa/2+sc,sa/2+sc,sa+sc+se),2,2) vcvm <- matrix(c(sa+sc+se,sa+sc,sa+sc,sa+sc+se),2,2) is it correct that these two differ? Also is the order of operations correct here? Maybe what you want is sa/(2+sc) (just a wild guess)? 3. Focus on accuracy first, then worry about efficiency ? but while we are at it that loop can definitely be eliminated, either by using apply (apply(D,1,FUN = muvcv.fun(mu = mu, vcv = vcv)) replace this by apply(D,1,FUN = muvcv.fun,mu = mu, vcv = vcv) as we discussed for passing arguments within apply) or by simply rewriting the whole expression in a matrix way (eg. X-mu is cbind(X[,1]-mean(X[,1], X[,2]-mean(X[, 2]), etc). 4. I don't see why the search space should be highly unregluar here (but you can check this which a little grid search), so optim should work. If not, try genoud(rgenoud), optimize() won't help you at all. 5. the error you are getting could be due to ill-conditioned matrices in vcvd. Check the condition numbers on the matrix, if it's ill- conditioned you cannot invert it using solve, but this should not happen unless your data is very highly collinear. hth, jens From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at lists.fas.harvard.edu] On Behalf Of Jeremy Hodgen Sent: Saturday, March 29, 2008 11:51 AM To: gov2001-l at lists.fas.harvard.edu Subject: [gov2001-l] Problems with Optim Dear All, Apologies for the long message, but we didn't think that our questions would make sense otherwise. We're having some problems getting optim to work with our model & we'd be grateful for advice. The paper we're trying to replicate involves a twin study and compares identical (MZ) and non-identical (DZ) twins in an attempt to estimate the effects of genetic and environmental effects. Sets of twins are modelled using a bivariate normal with identical and non-identical twins having different but related distributions - basically the model is that identical and non-identical twins share genetic causes 100% and 50% respectively. Everything else (environment) is assumed to be equivalent. The model is a simple one - it assumes 3 latent causal variables: genetic, shared environment and non-shared environment. We are (at least initially) interested in the effects of individual traits (e.g., Speaking at age 7, as measured by teacher assessment.) If we get past this, then we'll try and replicate their longitudinal modeling. We are interested in the proportions of these effects on the variance & covariance. So we've created a log-likelihood function based on this (see R code pasted below) and we want to produce a MLE for the three factors. But optim doesn't work - it produces the following error message: Error in solve.default(vcvm) : Lapack routine dgesv: system is exactly singular We've tried all the different methods that optim uses (except "L- BFGS-B") and all (except "SANN") produce the same error message (sometimes referring to vcvd.) We've left SANN going all night and it hasn't produced an answer. Now vcvd and vcvm are variance-covariance matrices and, our function does attempt to work out the inverse using solve, but the way in which we have set them up should not allow singular matrices to be produced (see the code.) However, we've tried a range of different starting values & we get the same errors. Our data has some odd features. After excluding some cases, there's roughly 2000 identical or MZ twins and roughly 3500 non-identical or DZ twins with lots of missing values. But many of the traits that we're interested in have a very small scale - so Speaking at age 7 is just a scale of 3 levels (1, 2 & 3) - and of course the trait for many of the twins (about 75% of them) are exactly the same. We can argue about whether this is a good model - but our initial job is to replicate the results of the paper. There are several other traits that are more straightforward (so 'g' or intelligence is measured as z-scores.) So the specific questions are: 1. Why is optim doing this and should we use something different? The help file on optim suggests that optim might be problematic for parameters with more than one dimension (which is really what we've got here since we use the 3 variance parameters to create a variance-covariance matrix)? The help files suggests using optimize, but it's not at all clear how we should use optimize - the arguments don't appear to us to be equivalent to optim (e.g., no par). Also we're uninterested in mu - should we actually just work with the actual mean for the trait from our data - which would presumably be faster. 2. Is there a way of avoiding the for loop in the code? It seems to us that this will be really inefficient & optim certainly takes a while to return the error message (> 2 minutes). We've tried using apply with a function we define ourselves: mu <- c(1,1) vcv <- matrix(c(1,0,0,1),2,2) D <- matrix(c(1,2,3,4),2,2) muvcv.fun <- function(x,mu,vcv) (x-mu) %*% vcv %*% (x-mu) apply(D,1,FUN = muvcv.fun(mu = mu, vcv = vcv)) But we get the following error message: Error in muvcv.fun(mu = mu, vcv = vcv) : argument "x" is missing, with no default So we're unsure how the syntax for a parameter should be for a function inside optim. We've tried various options that seem intuitive / sensible to us, although doubtless we're missing something obvious. 3. We have to then calculate the root mean square error of approximation (RMSEA). In the original paper, this is calculated using the software the authors use - Mx, which seems to be pretty exclusively used for twins studies, although RMSEA is fairly commonplace. The formula we've found for RMSEA is sqrt[(chi^2/df - 1)/(N-1)] but it's not very clear how to use this. We think we have 6 parameters (the variances and covariances for the identical and non-identical twins] - 3 for the three variance proportions that we're estimating, hence 3 df. Then N is the total sample size. But to get chi^2, do we take (obs-exp)^2/exp? Thanks Jeremy & Jill (Hohenstein) Here's the R-code: # Implement log-likelihood function for MZ & DZ twins - simple comparison of a single trait based ACE model & bivariate normal distribution ll.dbnmbn <- function(par,Yd,Ym,Xd,Xm){ sa <- exp(par[2]) # Re-parameterize R -> R+ as sigma^2 genetic (hereditary) variance sc <- exp(par[3]) # Shared environment variance se <- exp(par[4]) # Non-shared environment mu <- par[1] vcvd <- matrix(c(sa+sc+se,sa/2+sc,sa/2+sc,sa+sc+se),2,2) # varcov matrix for DZ twins vcvm <- matrix(c(sa+sc+se,sa+sc,sa+sc,sa+sc+se),2,2) # varcov matrix for MZ twins D <- rbind(Yd,Xd) # 2 x N(DZ) matrix of bivariate DZ twins trait scores M <- rbind(Ym,Xm) # 2 x N(MZ) matrix of bivariate MZ twins trait scores d <- 0 for (i in 1:length(Xd)){ d <- d + (D[,i]-mu)%*%solve(vcvd)%*%(D[,i]-mu) # sum((X-mu) x varcov^-1 x (X-mu)) for DZ twins } m <- 0 for (i in 1:length(Xm)){ m <- m + (M[,i]-mu)%*%solve(vcvm)%*%(M[,i]-mu) # sum((X-mu) x varcov^-1 x (X-mu)) for MZ twins } # Log-likelihood func based on bivariate normal - related but different for DZ and MZ twins out <- -length(Xd)*log(2*pi*sqrt(det(vcvd)))-0.5*d - length(Xm)*log(2*pi*sqrt(det(vcvm)))-0.5*m return(out) } fit.1 <- optim(par = c(0,2.3,1.5,13.5), ll.dbnmbn, Yd = Yd, Ym = Ym, Xd = Xd, Xm = Xm, method="BFGS", control=list(fnscale = -1), hessian = T) Dr Jeremy Hodgen Senior Lecturer in Mathematics Education King's College London Department of Education and Professional Studies Franklin-Wilkins Building Waterloo Bridge Wing 150 Stamford Street London SE1 9NH Tel: 020 7848 3102 Fax: 020 7848 3182 E-mail: jeremy.hodgen at kcl.ac.uk _______________________________________________ gov2001-l mailing list gov2001-l at lists.fas.harvard.edu http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l