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]]
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: rvaradhan at jhmi.edu ----- Original Message ----- From: "Doran, Harold" <HDoran at air.org> Date: Wednesday, September 29, 2010 10:36 am Subject: [R] nlminb and optim To: R-help <R-help at 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 at r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.