9
$\begingroup$

I have some time series data where the measured variable is discrete positive integers (counts). I want to test if there is an upward trend over time (or not). The independent variable (x) is in the range 0-500 and the dependent variable (y) is in the the range 0-8.

I thought that I answer this by fitting a regression of the form y = floor(a*x + b) using ordinary least squares (OLS).

How would I go about doing this using R (or Python)? Is there an existing package for it, or am I better off writing my own algorithm?

PS: I know this is not the ideal technique, but I need to do a relatively simple analysis that I can actually understand - my background is biology not maths. I know I am violating assumptions about error in measured variable, and independence of measurements over time.

$\endgroup$
  • 5
    $\begingroup$ Although it is mathematically natural to try a regression of this form, behind it lurks a statistical mistake: the error term will now be strongly correlated with the predicted value. That's a pretty strong violation of OLS assumptions. Instead, use a count-based technique as suggested by Greg Snow's reply. (I gladly upvoted this question, though, because it reflects some real thought and cleverness. Thank you for asking it here!) $\endgroup$ – whuber Jan 13 '12 at 21:37
11
$\begingroup$

You could fit the model you state using the nls (non-linear least squares) function in R, but as you said that will violate many of the assumptions and still probably will not make much sense (you are saying the predicted outcome is random around a step function, not integer values around a smoothly increasing relationship).

The more common way to fit count data is using Poisson regression using the glm function in R, the first example on the help page is a Poisson regression, though if you are not that familiar with statistics it would be best to consult with a statistician to make sure that you are doing things correctly.

If the value of 8 is an absolute maximum (impossible to ever see a higher count, not just that is what you saw) then you might consider proportional odds logistic regression, there are a couple of tools to do this in packages for R, but you really should get a statistician involved if you want to do this.

$\endgroup$
  • $\begingroup$ "you are saying the predicted outcome is random around a step function, not integer values around a smoothly increasing relationship" --- That is something I had not considered. In the end, I went with Poisson regression by glm. Its not the perfect choice, but "good enough" for what I needed. $\endgroup$ – afaulconbridge Jan 20 '12 at 14:40
10
$\begingroup$

$\def\lf{\lfloor}\def\rf{\rfloor}\def\pnorm{\mathrm{pnorm}}$It is plain that Greg’s suggestion is the first thing to try: Poisson regression is the natural model in many many concrete situations.

However the model you’re suggesting can occur for example when you observe rounded data: $$ Y_i = \lf ax_i + b + \epsilon_i \rf,$$ with iid normal errors $\epsilon_i$.

I think this is interesting to have a look on what can be done with it. I denote by $F$ the cdf of the standard normal variable. If $\epsilon \sim \mathcal N(0,\sigma^2)$, then $$\begin{align*} \mathbb P\left(\lf ax + b + \epsilon \rf = k\right) &= F\left({k-b+1-ax\over \sigma}\right) - F\left({k-b-ax\over \sigma}\right)\\ &= \pnorm(k+1-ax-b,sd=\sigma) - \pnorm(k-ax-b,sd=\sigma),\end{align*}$$ using familiar computer notations.

You observe data points $(x_i,y_i)$. The log likelihood is given by $$\ell(a,b,\sigma) = \sum_i \log\left( F\left({y_i-b+1-ax_i\over \sigma}\right) - F\left({y_i-b-ax_i\over \sigma}\right) \right).$$ This is not identical to least squares. You can try to maximize this with a numerical method. Here is an illustration in R:

log_lik <- function(a,b,s,x,y)
  sum(log(pnorm(y+1-a*x-b, sd=s) - pnorm(y-a*x-b, sd=s)));

x <- 0:20
y <- floor(x+3+rnorm(length(x), sd=3))
plot(x,y, pch=19)
optim(c(1,1,1), function(p) -log_lik(p[1], p[2], p[3], x, y)) -> r
abline(r$par[2], r$par[1], lty=2, col="red")
t <- seq(0,20,by=0.01)
lines(t, floor( r$par[1]*t+r$par[2]), col="green")

lm(y~x) -> r1
abline(r1, lty=2, col="blue");

rounded linear model

In red and blue, the lines $ax+b$ found by numerical maximization of this likelihood, and least squares, respectively. The green staircase is $\lf ax +b\rf$ for $a,b$ found from the max likelihood... this suggest that you could use least squares, up to a translation of $b$ by 0.5, and get roughly the same result; or, that least squares fit well the model $$ Y_i = [ a x_i + b +\epsilon_i], $$ where $[x] = \lf x + 0.5 \rf$ is the nearest integer. Rounded data are so often met that I am sure this is known and had been studied extensively...

$\endgroup$
  • 4
    $\begingroup$ +1 I love this technique and actually submitted a paper on it to a risk analysis journal a few years ago. (Some risk analysts are quite interested in interval-valued data.) It was rejected as being "too mathematical" for their audience. :-(. One tip: when using numerical methods, it's always a good idea to supply good starting values for the solution. Consider applying OLS to the raw data to obtain those values, then "polish" them with the numerical optimizer. $\endgroup$ – whuber Jan 14 '12 at 18:09
  • $\begingroup$ Yes, this is a good suggestion. In fact, in that case I choose remote values to emphasize that "it works", but in practice your suggestion would be the only solution to avoid to start from a very flat region, depending on the data... $\endgroup$ – Elvis Jan 14 '12 at 22:42

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.