Can you send your code and data as separate files so we can get it into R easily? ____________________________________________________________________
Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: "Doran, Harold" <hdo...@air.org> Date: Wednesday, September 29, 2010 10:36 am Subject: [R] nlminb and optim To: R-help <R-help@r-project.org> > I am using both nlminb and optim to get MLEs from a likelihood > function I have developed. AFAIK, the model I has not been previously > used in this way and so I am struggling a bit to unit test my code > since I don't have another data set to compare this kind of estimation > to. > > The likelihood I have is (in tex below) > > \begin{equation} > \label{eqn:marginal} > L(\beta) = \prod_{s=1}^N \int > \prod_{i=1}^K\frac{e^{x_{is}(\theta_s-\beta_i)}} > {x_{is}!e^{e^(\theta_s-\beta_i)}} f(\theta)d(\theta) > \end{equation} > > Where I view $\theta$ as a nuisance parameter and so I integrate it > out of the likelihood. The goal is to get parameter estimates for > $\beta$. The integral cannot be easily evaluated so I approximate it as: > > \begin{equation} > \label{eqn:marginal2} > L(\beta) \approx \prod_{s=1}^N \sum_{q=1}^Q > \prod_{i=1}^K\frac{e^{x_{is}(\theta_{q}-\beta_i)}} > {x_{is}!e^{e^(\theta_{q}-\beta_i)}} w_q > \end{equation} > > \noindent where $\theta_{q}$ is the node at quadrature point $q = 1, > \ldots, Q$ and $w_q$ is the weight at quadrature point $q$. > > For now, I am assuming $f(\theta)$ is Uniform but this may change and > that is not a major issue. Now, I have written a function using both > nlminb and optim. I have copied my function below where the arguments > are the data set, Q for the number of quadrature points, and an option > for providing starting values. I am running into some flow control > problems and so I hack it a bit and fix a lower bound as you may see > in the function. > > Here is the issue at hand. First, nlminb and optim give different > parameter estimates. Optim tells me it converged but nlminb tells me I > get relative convergence. If I use different number of quadrature > points in optim, I get very different parameter estimates and this > should not happen. But, varying the number of quadrature points in > nlminb yields the same parameter estimates, but I am not sure the > model converged. > > I have also provided the data set I am working with below as well as > sessionInfo(). There may be many issues going on here and am grateful > for any reactions you all may have. I *think* my likelihood is written > properly and I *think* I am handling flow control issues in a > relatively decent way. I may be misinterpreting optim or nlminb reports. > > In any event, should anyone take the time to review this and offers > pointers I would be grateful. > Harold > > > > fit <- function(data, Q, startVal = NULL, ...){ > if(!is.null(startVal) && length(startVal) != > ncol(data) ){ > stop("Length of argument startVal not > equal to the number of parameters estimated") > } > if(!is.null(startVal)){ > startVal <- startVal > } else { > startVal <- rep(0, ncol(data)) > } > #qq <- gauss.quad.prob(Q, dist = 'normal', mu = 0, sigma=1) > qq <- gauss.quad.prob(Q, dist = 'uniform', alpha = -5, > beta=5) > rr1 <- matrix(0, nrow = Q, ncol = nrow(data)) > data <- as.matrix(data) > L <- nrow(data) > C <- ncol(data) > fn <- function(b){ > b <- b[1:C] > for(j in 1:Q){ > for(i in 1:L){ > > rr1[j,i] <- prod((exp(data[i,]*(qq$nodes[j]-b))) / > (factorial(data[i,]) * exp(exp(qq$nodes[j]-b)))) * > qq$weights[j] > } > } > rr1[rr1==0] <- 1e-10 > -sum(log(colSums(rr1))) > } > #opt <- optim(startVal, fn, method = "BFGS", hessian = > FALSE) > opt <- nlminb(startVal, fn) > #list("coefficients" = opt$par, "LogLik" = -opt$value, > "Std.Error" = sqrt(diag(solve(opt$hessian)))) > #list("coefficients" = opt$par, "LogLik" = -opt$value) > opt > } > > fit(data, Q = 10) > > > > > data > cindyR cjR ohsR sdhpR > [1,] 24 14 6 9 > [2,] 19 14 6 10 > [3,] 19 14 6 10 > [4,] 17 18 6 8 > [5,] 19 17 5 8 > [6,] 17 17 6 9 > [7,] 13 15 5 11 > [8,] 19 13 6 8 > [9,] 21 10 4 11 > [10,] 16 15 5 9 > [11,] 14 16 5 10 > [12,] 14 13 5 10 > [13,] 17 15 5 8 > [14,] 16 18 4 8 > [15,] 16 14 5 8 > [16,] 18 13 5 8 > [17,] 15 14 6 7 > [18,] 16 17 5 7 > [19,] 18 12 5 8 > [20,] 18 9 4 9 > [21,] 17 9 4 10 > [22,] 18 12 4 8 > [23,] 15 15 5 7 > [24,] 16 9 5 8 > [25,] 18 11 5 7 > [26,] 18 10 4 9 > [27,] 15 15 4 7 > [28,] 16 11 4 9 > [29,] 11 16 5 8 > [30,] 11 16 5 8 > [31,] 16 13 4 8 > [32,] 13 15 5 7 > [33,] 17 17 2 9 > [34,] 18 8 3 9 > [35,] 13 14 4 8 > [36,] 16 11 4 8 > [37,] 16 11 4 8 > [38,] 15 8 5 9 > [39,] 14 15 4 8 > [40,] 15 14 4 7 > [41,] 14 11 4 8 > [42,] 13 12 4 9 > [43,] 15 13 3 8 > [44,] 12 17 3 8 > [45,] 20 0 4 8 > [46,] 14 10 4 9 > [47,] 11 15 5 7 > [48,] 14 12 4 8 > [49,] 11 14 4 8 > [50,] 13 11 5 7 > [51,] 12 13 5 6 > [52,] 16 7 5 7 > [53,] 10 13 4 8 > [54,] 15 12 3 7 > [55,] 16 7 4 8 > [56,] 16 11 3 7 > [57,] 12 15 4 7 > [58,] 19 2 4 8 > [59,] 13 13 4 7 > [60,] 15 12 3 7 > [61,] 14 15 3 7 > [62,] 19 10 2 8 > [63,] 13 7 5 7 > [64,] 16 7 4 7 > [65,] 16 2 5 7 > [66,] 12 17 3 7 > [67,] 13 6 4 9 > [68,] 15 7 4 7 > [69,] 13 11 3 8 > [70,] 15 7 4 7 > [71,] 12 14 4 6 > [72,] 14 6 4 8 > [73,] 15 6 4 7 > [74,] 13 11 4 7 > [75,] 15 8 4 7 > [76,] 13 1 6 7 > [77,] 12 13 3 7 > [78,] 13 6 4 8 > [79,] 16 10 2 7 > [80,] 11 20 3 6 > [81,] 12 9 4 7 > [82,] 16 6 3 8 > [83,] 18 6 2 8 > [84,] 12 12 4 6 > [85,] 14 10 2 7 > [86,] 15 1 4 8 > [87,] 12 9 3 8 > [88,] 12 10 4 6 > [89,] 14 6 3 7 > [90,] 15 0 5 7 > [91,] 17 3 1 9 > [92,] 14 7 2 7 > [93,] 13 10 2 8 > [94,] 13 5 4 7 > [95,] 18 2 2 7 > [96,] 15 5 3 7 > [97,] 11 3 5 7 > [98,] 11 6 4 6 > [99,] 16 1 3 8 > [100,] 13 7 3 6 > [101,] 13 4 4 6 > [102,] 15 9 2 6 > [103,] 13 4 3 7 > [104,] 14 4 3 7 > [105,] 17 0 2 8 > [106,] 13 3 2 8 > [107,] 9 7 4 7 > [108,] 11 4 4 6 > [109,] 7 12 3 6 > [110,] 9 7 4 6 > [111,] 14 2 1 9 > [112,] 13 15 0 6 > [113,] 9 5 3 8 > [114,] 11 5 4 6 > [115,] 13 2 2 8 > [116,] 9 4 4 6 > [117,] 14 2 2 6 > [118,] 13 8 0 8 > [119,] 11 7 3 5 > [120,] 16 1 1 7 > [121,] 13 1 3 6 > [122,] 11 4 2 7 > [123,] 11 1 2 8 > [124,] 13 2 2 6 > [125,] 11 5 2 6 > [126,] 14 0 2 7 > [127,] 11 7 2 6 > [128,] 8 8 3 6 > [129,] 10 4 3 6 > [130,] 13 0 2 6 > [131,] 10 2 2 7 > [132,] 10 1 3 6 > [133,] 11 2 2 6 > [134,] 11 0 2 6 > [135,] 10 1 3 6 > [136,] 11 1 2 6 > [137,] 14 1 2 5 > [138,] 10 1 1 7 > [139,] 12 0 2 6 > [140,] 10 0 2 6 > [141,] 11 0 0 6 > > > > sessionInfo() > R version 2.10.1 (2009-12-14) > i386-pc-mingw32 > > locale: > [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United > States.1252 LC_MONETARY=English_United States.1252 > [4] LC_NUMERIC=C LC_TIME=English_United States.1252 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] MiscPsycho_1.6 statmod_1.4.6 lattice_0.17-26 gdata_2.8.0 > > loaded via a namespace (and not attached): > [1] grid_2.10.1 gtools_2.6.2 tools_2.10.1 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.