hypothesis testingnonlinear regressionr
I have a set of results that take the following form:
σ Standard error
0,06 0,50
0,91 >10
0,04 0,11
0,84 0,56
0,44 0,10
0,03 0,09
0,30 >10
0,19 0,07
0,04 1,15
Where σ is a coefficient in a non-linear model, which I've estimated directly using an iterative algorithm (BFGS).
If I want to test a hypothesis of the type σ=0, which test would be appropriate? We haven't covered non-linear optimization methods in my classes, and I'm really at loss.
The weights should equal the counts, because those will be inversely proportional to the variances of the errors. Specifically, the model for the data $(x_i, y_i, n_i)$ is
$$y_i \sim \lambda \Phi((\log(x_i) - \mu)/\sigma + \varepsilon_i$$
with $\mu, \sigma \gt 0,$ and $\lambda \gt 0$ the parameters and $\varepsilon_i$ are independent random variables with zero means and variances
$$\text{Var}(\varepsilon(i)) = \sigma^2 / n_i$$
where $n_i$ are the counts.
The fit to the logarithm of $x$ is visually ok:
In this figure the x-axis is on a logarithmic scale, the point symbols have areas proportional to the counts (so that large circles will have more influence in the fitting than small ones), and the red line is a least-squares fit. It is clear the model is not really appropriate: the residuals for smaller values of $y$ tend to be small, regardless of the counts. Possibly the sum of squares of relative errors should be minimized to obtain an appropriate fit.
It is evident that the fit is poor for the largest $x$, but those also have small counts.
The R code with (my version of) the data and the fitting and plotting procedures follows.
y <- c(1, 1, 2, 1, 2, 1, 3, 4, 22, 30, 44, 58, 68, 69,
71, 72, 75, 72, 80, 78, 87, 86, 80, 82, 92, 90, 85, 61, 38, 36) / 100
x <- ceiling(exp(seq(log(20), log(500), length.out=length(y))))
counts <- c( 10, 3, 17, 20, 38, 31, 44, 55, 58, 68, 77,
82, 86, 82, 77, 75, 70, 65, 68, 51, 47, 41, 38, 30, 22, 14, 9, 4, 2, 1)
#
# The least-squares criterion.
# theta[1] is a location, theta[2] an x-scale, and theta[3] a y-scale.
#
f <- function(theta, x=x, y=y, n=counts)
sum(n * (y - pnorm(x, theta[1], theta[2]) * theta[3])^2) / sum(n)
#
# Perform a count-weighted least-squares fit.
#
xi = log(x)
fit <- optim(c(median(xi), sd(xi), max(y) * sd(xi)), f, x=xi, y=y, n=counts)
#
# Plot the result.
#
par(mfrow=c(1,1))
plot(x, y, log="x", xlog=TRUE, pch=19, col="Gray", cex=sqrt(counts/12))
points(x, y, cex=sqrt(counts/10))
curve(fit$par[3] * pnorm(log(x), fit$par[1], fit$par[2]),
from=10, to=1000, col="Red", add=TRUE)
I fit your data with AD Model Builder. The R code has the right std devs (almost) but poor parameter estimates. The C code has good parameter estimates but the wrong std devs You need to multiply by sqrt(9/7) to convert from fisher information to nls stdevs
index name value std.dev
1 N 4.7079e+00 1.7120e-01
2 rate 4.7979e-01 1.9154e-02
Best Answer
Hypothesis testing and confidence intervals on non-linear regression models are based on asymptotic theory. To construct a test, you may extract the diagonal element of the covariance matrix that corresponds to your estimate and as usual form the t-ratio. You can then use the critical values of the standard normal distribution to reach a decision.
One warning is in order though. As these tests are based on large sample theory, it's not prudent to pay much attention to them when your sample is quite small. There is no clear guideline as to "how large is sufficiently large" but personally if it's less than 50-60 observations, I would be very suspicious.
Hope this helps.