Jeremy,
you need to pass the other function arguments to genoud (just as in optim), otherwise it
doesn?t know where to look. Try:
genoud(ll.dbnmbn, nvars = 3, max =T, hessian = T , mu=mu,D=D,M=M)
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: Sunday, March 30, 2008 7:21 PM
To: gov2001-l at
lists.fas.harvard.edu
Subject: Re: [gov2001-l] Problems with Optim
Jens,
Thanks. We're trying to fit two highly inter-related but not identical bivariate
normals to our data - see attached. Although maybe we're thinking about it
incorrectly. There's some info on the site our paper comes from to suggest that there
may be several local maxima - and this process requires patience!
We can't get genoud to work though - & get this error message with the following
code:
Maximization Problem.
Error in apply(M, 1, FUN = fnxvx, y = mu, S = vcvm) :
argument "M" is missing, with no default
R-code:
D <- Db # Toy data
M <- Mb
mu <- DMmean.b
genoud(ll.dbnmbn, nvars = 3, max =T, hessian = T)
And the function is:
# Implement log-likelihood function for MZ & DZ twins - simple comparison of a single
trait based ACE model & bivariate normal distribution
fnxvx <- function(X,y,S) (X-y)%*%solve(S)%*%(X-y) # Function to calculate (X-mu) x
var-cov^-1 x (X-mu)
# X <- c(11,7); y <-c(1,2); S <- matrix(c(1,0,0,1),2,2); fnxvx(X,y,S) # Test
fnxvx
ll.dbnmbn <- function(par,mu,D,M){
sa <- exp(par[1]) # Re-parameterize R -> R+ as sigma^2 genetic (hereditary)
variance (sigma^2[A])
sc <- exp(par[2]) # Shared environment variance (sigma^2[C])
se <- exp(par[2]) # Non-shared environment variance (sigma^2[E])
vcvm <- matrix(c(sa+sc+se,sa+sc,sa+sc,sa+sc+se),2,2) # varcov matrix for MZ
(identical) twins
vcvd <- matrix(c(sa+sc+se,sa/2+sc,sa/2+sc,sa+sc+se),2,2) # varcov matrix for DZ
(non-identical) twins
m <- sum(apply(M,1,FUN=fnxvx,y=mu,S=vcvm)) # sum((X-mu) x varcov^-1 x (X-mu)) for MZ
(identical) twins
d <- sum(apply(D,1,FUN=fnxvx,y=mu,S=vcvd)) # sum((X-mu) x varcov^-1 x (X-mu)) for DZ
(non-identical) twins
# Log-likelihood func based on bivariate normal - related but different for DZ and MZ
twins
out <- -nrow(M)*log(det(vcvm))-m - nrow(D)*log(det(vcvd))-d
return(out)
}
And our toy data:
# Create some alternative toy data based on mu = 100, sa = 6, sc = 2, se = 2
set.seed(12345)
Db <- mvrnorm(1000,c(100,100),matrix(c(10,8,8,10),2,2)); head(Db); nrow(Db)
Mb <- mvrnorm(1000,c(100,100),matrix(c(10,7,7,10),2,2)); head(Mb); nrow(Mb)
DMmean.b <- (apply(Db,2,mean) + apply(Mb,2,mean))/2; DMmean
ll.dbnmbn(log(c(6,2,2)),100,Da,Ma) # Does the function still work? Yes
true.b <- c(log(6), log(2), log(2)); true.b # Numbers optim SHOULD return
With optim we get different but pretty close results from different starting values - but
some starting values produce the solve/singular error we got earlier.
Jeremy