# I’ll take my NLS with weights, please…

January 13, 2014
By

(This article was first published on Rmazing, and kindly contributed to R-bloggers)

Today I want to advocate weighted nonlinear regression. Why so?
Minimum-variance estimation of the adjustable parameters in linear and non-linear least squares requires that the data be weighted inversely as their variances $w_i \propto \sigma^{-2}$. Only then $\hat{\beta}$ is the BLUE (Best Linear Unbiased Estimator) for linear regression and nonlinear regression with small errors (http://en.wikipedia.org/wiki/Weighted_least_squares#Weighted_least_squares), an important fact frequently neglected, especially in scenarios with heteroscedasticity.
The variance of a fit $s^2$ is also characterized by the statistic $\chi^2$ defined as followed:
$\chi^2 \equiv \sum_{i=1}^n \frac{(y_i - f(x_i))^2}{\sigma_i^2}$
The relationship between $s^2$ and $\chi^2$ can be seen most easily by comparison with the reduced $\chi^2$: $\chi_\nu^2 = \frac{\chi^2}{\nu} = \frac{s^2}{\langle \sigma_i^2 \rangle}$
whereas $\nu$ = degrees of freedom (N – p), and $\langle \sigma_i^2 \rangle$ is the weighted average of the individual variances. If the fitting function is a good approximation to the parent function, the value of the reduced chi-square should be approximately unity, $\chi_\nu^2 = 1$. If the fitting function is not appropriate for describing the data, the deviations will be larger and the estimated variance will be too large, yielding a value greater than 1. A value less than 1 can be a consequence of the fact that there exists an uncertainty in the determination of $s^2$, and the observed values of $\chi_\nu^2$ will fluctuate from experiment to experiment. To assign significance to the $\chi^2$ value, one can use the integral probability $P_\chi(\chi^2;\nu) = \int_{\chi^2}^\infty P_\chi(x^2, \nu)dx^2$ which describes the probability that a random set of $n$ data points sampled from the parent distribution would yield a value of $\chi^2$ equal to or greater than the calculated one. This can be calculated by 1 - pchisq(chi^2, nu) in R.

To see that this actually works, we can Monte Carlo simulate some heteroscedastic data with defined variance as a function of $y$-magnitude and compare unweighted and weighted NLS.
First we take the example from the documentation to nls and fit an enzyme kinetic model:
 DNase1 <- subset(DNase, Run == 1) fm3DNase1 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1)) 
Then we take the fitted values $\hat{y}$ (which are duplicated because of the initial replicates), create a new unique dataset on which we create 20 response values $y_i$ for each concentration $x$ sampled from a normal distribution with 2% random heteroscedastic gaussian noise as a function of the value’s magnitude $y_i = \hat{y} + \mathcal{N}(0, 0.02 \cdot \hat{y})$:
 FITTED <- unique(fitted(fm3DNase1)) DAT <- sapply(FITTED, function(x) rnorm(20, mean = x, sd = 0.02 * x)) matplot(t(DAT), type = "p", pch = 16, lty = 1, col = 1) lines(FITTED, col = 2) 

Now we create the new dataframe to be fitted. For this we have to stack the unique $x$– and $y_i$-values into a 2-column dataframe:
 CONC <- unique(DNase1$conc) fitDAT <- data.frame(conc = rep(CONC, each = 20), density = matrix(DAT))  First we create the unweighted fit:  FIT1 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = fitDAT, start = list(Asym = 3, xmid = 0, scal = 1))  Then we fit the data with weights $w = 1/var(y)$. IMPORTANT: we need to replicate the weight values by 20 in order to match the data length.  VAR <- tapply(fitDAT$density, fitDAT$conc, var) VAR <- rep(VAR, each = 20) FIT2 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = fitDAT, weights = 1/VAR, start = list(Asym = 3, xmid = 0, scal = 1))  For calculation of $\chi^2_\nu$ and its corresponding p-value, we use the fitchisq function of my ‘qpcR’ package:  library(qpcR) > fitchisq(FIT1)$chi2 [1] 191.7566 $chi2.red [1] 1.22138$p.value [1] 0.03074883

 

> fitchisq(FIT2) $chi2 [1] 156.7153$chi2.red [1] 0.9981866 \$p.value [1] 0.4913983 

Now we see the benefit of weighted fitting: Only the weighted model shows us with it’s reduced chi-square value of almost exactly 1 and its high p-value that our fitted model approximates the parent model. And of course it does, because we simulated our data from it…

Cheers,
Andrej

Filed under: nonlinear regression Tagged: chi-square, nonlinear, regression, weights

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...