Which results should I trust between command “Hessian” and “numericHessian”?

You have not given a reproducible example, but I'll try anyway.

library(bbmle)
 x <- 0:10
 y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
 d <- data.frame(x,y)
 LL <- function(ymax=15, xhalf=6)
     -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
 fit <- mle2(LL)
 cc <- coef(fit)

Here are the finite-difference estimates of the Hessians (matrices of second derivatives) of the negative log-likelihood function at the MLE: inverting these matrices gives an estimate of the variance-covariance matrices of the parameters.

 library(numDeriv)
 hessian(LL,cc)
 ##               [,1]          [,2]
 ## [1,]  1.296717e-01 -1.185789e-15
 ## [2,] -1.185789e-15  4.922087e+00

 library(maxLik)
 numericHessian(LL, t0=cc)
 ##           [,1]     [,2]
 ## [1,] 0.1278977 0.000000
 ## [2,] 0.0000000 4.916956

So for this relatively trivial example, numDeriv::hessian and maxLik::numericHessian give very similar results. So there must be something you haven't shown us, or something special about the numerics of your problem. In order to proceed further, we need a reproducible example please ...

dat <- matrix(c(137,201,122,342,142,111,
                171,126,134,123,823,876,
                634,135,541,214,423,142),
        byrow=TRUE,ncol=2)
f <- function(x){
    w1 <- tail(x/(rowSums(x)),1)
    sqrt(w1%*%cov(cbind(x))%*%t(w1))
}
p <- t(dat[1,1:2,drop=FALSE])
f(p)  ## 45.25483
numDeriv::hessian(f,p)
##              [,1]          [,2]
## [1,] -3.217880e-15 -1.957243e-16
## [2,] -1.957243e-16  1.334057e-16
maxLik::numericHessian(f,t0=p)
##            [,1]        [,2]
## [1,] 0.007105427 0.007105427
## [2,] 0.007105427 0.000000000

OK, these clearly disagree. I'm not sure why, but in this particular case we can analyze what you're doing and see which one is right:

• since your input matrix is a single column, x/rowSums(x) is a vector of ones, so the last element (w1 <- tail(...,1)) is just 1.
• so your expression reduces to sqrt(cov(cbind(x))). Again, since x is a one-column matrix, cov() is just the variance, and sqrt(cov(.)) is just the standard deviation, or the norm of the vector.
• the variance is a quadratic function of any element's deviation from the mean, and so the standard deviation is more or less linear in the deviation from the mean (except at zero), so we would expect the second derivatives to be zero. So it looks like numDeriv::hessian is giving the right answer

We can also confirm this by increasing eps for numericHessian:

maxLik::numericHessian(f,t0=p,eps=1e-3)
##      [,1]          [,2]
## [1,]    0  0.000000e+00
## [2,]    0 -7.105427e-09

The bottom line is that numDeriv uses a more accurate (but slower) method, but you can get reasonable answers from numericHessian if you're careful.