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]]

