Standard errors for non-linear least squares in R

Question

I have a question on standard errors for non-linear least squares in R. With the built-in function NLS and a hand-made function I get different SE and I don't understand why. I will try to expose below a simpler and replicable example of my original problem.

I want to estimate the mean $\mu$ and the standard deviation $\sigma$ of a normal cumulative distribution function.

Let's simulate the data for this example in R:

set.seed(12)
mu <- 5
sigma <- 2
x <- runif(100, min = 0, max = 10)
y <- pnorm(x + rnorm(1), mu, sigma)
df <- data.frame(y,var1)

I assume a true value for $\mu$ of 5 and 2 for $\sigma$. I add a normal noise to y which is a cdf of x.

I first fit a non-linear least squares with the built-in function NLS:

res1 <- nls(y ~ pnorm(x, mean = mu, sd = sigma), 
            data = df, start = c(mu = 1,sigma = 1), algorithm = "port")
summary(res1)

I get the following result that I was expected

Formula: y ~ pnorm(x, mean = mu, sd = sigma)

Parameters:
                    Estimate             Std. Error            t value
mu    5.04268491217399539295 0.00000000000000003347 150682292747528064
sigma 2.00000000000000000000 0.00000000000000004567  43795073347627272
                 Pr(>|t|)    
mu    <0.0000000000000002 ***
sigma <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.00000000000000003995 on 98 degrees of freedom

Algorithm "port", convergence message: X-convergence (3)

Secondly, I tried my own LS estimation

f <- function(par){
  mu = par[1]
  sigma = par[2]
  temp <- pnorm(df$x, mean = mu, sd = sigma)
  sum1 <- sum((df$y-temp)^2)
  return(sum1)
}

res2 <- optim(c(1,1),f,method = "L-BFGS-B",
              lower=c(0,0), upper = c(100,100), hessian = TRUE)
print(res2)

I get similar result for the coefficients:

$par
[1] 5.042685 1.999999

$value
[1] 0.000000000001761186

$counts
function gradient 
      20       20 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$hessian
           [,1]       [,2]
[1,]  2.8835336 -0.2273554
[2,] -0.2273554  1.5485147

But when I compute the Standard Errors, I don't get the same result with the NLS:

sqrt(diag(res2$hessian))
[1] 1.698097 1.244393

Does someone know why ?

Thank you !

Alexis

Answer 1

I am almost certain that you are running into the same problem as here.

The thing is that optim is agnostic of any kind of least squares modeling, so you just get the pure hessian out. The std error has been scaled in nls.