How to fit a regression like $y=\lfloor ax+b \rfloor$ in R? 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.
 A: 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. 
A: $\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");


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... 
