I want to test how well my data fit the particular power law:


where $b$ , for physical reasons, should equal exactly $-0.5$. I would like to find the probability that the data do not obey this, or more correctly what the chance is of obtaining the data if it did equal $-0.5$. Here's some typical data with estimated $2\sigma$ uncertainties:

   x =(   5,  15,  13,  60, 125, 505)
delx =(   3,  14,   8,   2,   2,   3)
   y =( 415, 316, 168,  59,  29,   5)
dely =(  33,   3,   6,   4,   3,   1)

I hate to fit a straight line on a loglog plot in this neighbourhood, but the errors are so poorly estimated it hardly matters. The logging also biases the calculated power, but that should be ok as it just means that the $p$-value will be a maximum. The problem is that most of the error is in the $x$-axis-- is it permissible to swap $x$ and $y$ so a weighted linear regression takes account of the bulk of the error? I.e. instead of

$\log(y)=\log(a) + b \cdot log(x)$

use the model

$\log(x)=\frac{1}{b} \log(y) -\frac{1}{b} \log(a)$

Does a $t$-test on the power and its resulting $p$-value actually mean anything when there are so many violated assumptions? Also, would leave-one-out-cross-validation be a defence against the small sample?

  • 1
    $\begingroup$ The most critical aspect seems to be the errors in $x$ here. If, as you suggest, the errors in $y$ are comparatively very small you could fit the model $\text{E}(x) = (y/a)^{1/b} = \alpha y^\beta$ (where $\alpha = a^{-1/b}$ and $\beta = 1/b$) - that is, simply fit a power law with the roles of x and y swapped, and then back out parameter estimates (you should carefully consider your model for the variance of $x$ though; constant variance may not be suitable on any of the scales mentioned). If both variables have substantive error, you may need to consider errors-in-variables type models. $\endgroup$ – Glen_b Nov 12 '13 at 23:26

Because there are errors in both variables, we do not know the true values of the $(x,y)$ data associated with each observation. Let us suppose that each of $n=6$ independent observations results from measurement errors independently made in the two coordinates $(\xi, \eta) = (\xi, f(\xi;\alpha,\beta))$ with $f(\xi; \alpha, \beta) = \alpha \xi^\beta$. When those errors have Gaussian distributions and their standard deviations are proportional to the stated uncertainties $\sigma_i$ (in the first coordinate) and $\tau_i$ (in the second coordinate), the negative log likelihood of the data $(x_i, y_i)$ is

$$-\log(\Lambda(\alpha, \beta, (\xi_i))) = C+\frac{1}{2\lambda^2} \sum_i^n \left(\left(\frac{x_i - \xi_i}{\sigma_i}\right)^2 + \left(\frac{y_i - f(\xi_i; \alpha, \beta)}{\tau_i}\right)^2\right)$$

where the $\xi$ are the correct first coordinates and $C = n\log(2\pi) - 2 n \log(\lambda) + \sum_i^n \log(\sigma_i\tau_i)$ does not depend on the $n+2$ parameters $\alpha, \beta,$ and $(\xi_i)$.

The maximum likelihood estimates of $\alpha, \beta, (\xi_i)$ for these data (using $\lambda=1/2$) are $\hat\alpha = 7.77197$, $\hat\beta=-0.92872$, and $(\xi_i) = (6.5, 8.8, 17.4, 58.8, 124.9, 505.4)$, whence $(f(\xi_i; \alpha, \beta))$ = $(417.4, 316.0, 167.5, 54.0, 26.8, 7.3)$. To test whether $\beta=-1/2$ is consistent with the data, maximize the log likelihood with respect to that constraint. twice the difference in values at the optima under the null hypothesis has a $\chi^2(1)$ distribution. For these data it equals $1035,$ which is enormous: clearly $\beta\ne -1/2$.


In this figure, line segments connect the data to their fitted values for both fits. Symbol for the data are scaled in inverse proportion to the geometric means of the $\sigma_i$ and $\tau_i$, so that the more precise data are shown with larger circles. The slope of the data on this log-log plot is almost $-1$, which is far from the hypothesized value of $-1/2$. The best fit with a slope of $-1/2$ is shown in red. Evidently it posits tremendously large measurement errors occurred in both coordinates, far larger than suggested by the $(\sigma_i)$ and $(\tau_i)$.

If you mistrust the stated uncertainties, compensate by increasing $\lambda$. (In the R code below, include an option of the form , lambda=2 in the calls to optim.) Optimal values of $\lambda$ could be found via cross-validation to get a reasonable match between the residuals and the scaled standard deviations $\sigma_i\hat\lambda$ and $\tau_i\hat\lambda$, but I have not carried that out because the conclusion is so strong anyway. (Moreover, changing $\lambda$ will not change either of the fits.)

R code to perform the calculations and display the figure:

x.data <- matrix(c(5,15,13,60,125,505,
              3,14,8,2,2,3), ncol=2)
y.data <- matrix(c(415,316,168,59,29,5,
              33,3,6,4,3,1), ncol=2)
# Functional form of the fit.
f <- function(abt) {
  x <- abt[-(1:2)]
  exp(abt[2] * log(x) + abt[1])
# Negative log likelihood.
ll <- function(abt, x=x.data, y=y.data, lambda=1/2) {
  a <- abt[1]; b <- abt[2]; x.0 <- abt[-(1:2)]
  y.0 <- f(abt)
  sum((((x[ ,1] - x.0)/x[ ,2])^2 + ((y[ ,1] - y.0)/y[ ,2])^2) / (2*lambda^2) + 
        log(x[, 2]*y[, 2]/(lambda^2)))
# Constrained log likelihood.
ll.0 <- function(abt, slope, ...) {
  abt[2] <- slope
  ll(abt, ...)
# Brute force fitting using ML.
ab.0 <- coef(lm(log(y.data[ ,1]) ~ log(x.data[ ,1]))) # Initial LS estimates
abt.0 <- c(ab.0, x.data[ ,1])
fit <- optim(abt.0, ll, control=list(maxit=20000))
fit.0 <- optim(fit$par, ll.0, control=list(maxit=20000), slope=-1/2)
fit.0$par[2] <- -1/2
# Plot of data and the fits.
plot(x[,1], y[,1], xlim=c(1/2,600), ylim=c(2,600), log="xy",
     cex = 4 / sqrt(x[,2]*y[,2]),
     main="Data and Fits", xlab="X", ylab="Y")
segments(x[,1], y[,1], fit.0$par[-(1:2)],  f(fit.0$par), lty=3, col="#e0404080")
segments(x[,1], y[,1], fit$par[-(1:2)],  f(fit$par), lty=3, col="#4040e080")
points(x[,1], y[,1], xlim=c(1/2,600), ylim=c(2,600),
     cex = 4 / sqrt(x[,2]*y[,2]), pch=19, col="#80808080")
points(fit$par[-(1:2)], f(fit$par), pch=19, col="Blue")
points(fit.0$par[-(1:2)],  f(fit.0$par), pch=3, col="Red")
abline(c(1/log(10), 1)*coef(lm(log(f(fit$par)) ~ log(fit$par[-(1:2)]))), col="Blue")
abline(c(1/log(10), 1)*coef(lm(log(f(fit.0$par)) ~ log(fit.0$par[-(1:2)]))), col="Red")
legend(x="bottomleft", legend=c("Data", "Fit", "Fit (slope = -1/2)"),
       col=c("Gray", "Blue", "Red"),
       pch=c(19, 1, 3))
# Chi-squared statistic to test beta=-1/2.
(chisq <- 2*(fit.0$value - fit$value))
pchisq(chisq, 1, lower.tail=FALSE)
  • $\begingroup$ Thanks. Is there a standard reference textbook on this method, for more background? $\endgroup$ – reftt Jan 13 '14 at 11:05
  • $\begingroup$ If by "this method" you mean Maximum Likelihood, the answer is that there are loads of books: you can take your pick from elementary texts through advanced research. If you mean this particular application of ML to a nonlinear weighted errors-in-variables regression, I haven't seen exactly this approach, but it's simple and natural enough I'm sure others have used it before. $\endgroup$ – whuber Jan 13 '14 at 15:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.