# How to perform a non-linear regression pixel by pixel?

Using nlsLM from Package ‘minpack.lm’ is straightforward for a simple example like this:

``````MODEL:y=(exp(a*x+b*z+c)+d)^f
x=c(0.5,0.3,0.2,0.4)
z=c(0.1,0.6,1,0.9)
y=c(0.2,0.3,0.9,0.9)
fit=nlsLM(y~(exp(a*x+b*z+c)+d)^f,start = list(a = -.03, b = 0.5, c = 1,d=0.02,f=0.003))
``````

And no problem with that. My real data are multiple rasters for one year.So I have one year of data for x,y,z but as rasters that contain pixels,I want to do a non-linear regression pixel by pixel(each regression for one pixel is not related to the next pixel). example data:

``````   r <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
r1 <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
r2 <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))
st1=stack(r,r1,r2)

rl <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
rS <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
rT <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))
st2=stack(rl,rS,rT)

re <- raster(nrows=10, ncols=10); r <- setValues(r, 1:ncell(r))
ru <- raster(nrows=10, ncols=10);r1 <- setValues(r1, 1:ncell(r))
rg <- raster(nrows=10, ncols=10);r2 <- setValues(r1, 1:ncell(r))
st3=stack(re,ru,rg)
``````

so first x would be the first pixel of st1(for the whole temporal period)

so first y would be the first pixel of st2(for the whole temporal period)

so first z would be the first pixel of st3(for the whole temporal period)

do regression for this pixel then do the same for all other pixels

As a result I get a map (raster) 10*10 for each best parameter(a,b,c,d,f)

I hope it is clear and you have any help!

• There may be some underlying problems to resolve before doing any coding. In particular, you are trying to fit a five-parameter model to just four data points (or only three, in the case of the rasters!). Another issue--which won't go away even with more data--is that your model is not identifiable, because `exp(a*x+b*z+c)^d` is algebraically the same as `exp(d*a*x+d*b*z+d*c)`, which means only four of the five parameters can ever be determined. – whuber Mar 10 '15 at 20:45
• Would the regression example at the bottom of this page help? – MickyT Mar 10 '15 at 22:27
• @whuber Thanks for your comment.The model is not important for the moment but the principle of doing non-linear regression for multiple rasters is the most important. So you can please use whatever model to show we can do this? – usersam Mar 11 '15 at 8:21

This should do what you want.

The model I have used is not the one you posted and it probably doesn't make any sense at all, however it does demonstrate the principle.

``````# Set up the rasters
r1 <- r2 <- r3 <- r4 <- r5 <- r6 <- raster(nrows=10, ncols=10);
# Populate them with some values
r1 <- setValues(r1,runif(100,min=1,max=100));
r2 <- setValues(r2,runif(100,min=1,max=100));
r3 <- setValues(r3,runif(100,min=1,max=100));
r4 <- setValues(r4,runif(100,min=1,max=100));
r5 <- setValues(r5,runif(100,min=1,max=100));
r6 <- setValues(r6,runif(100,min=1,max=100));
# Stack them
st1 <- stack(r1,r2,r3);
st2 <- stack(r4,r5,r6);
# Set up the function
test <- function(r) {
x <- r[1:3];
y <- r[4:6];
result <- c(NA,NA);
try(result <- c(coef(nlsLM(y ~ a + b * x, start = list(a = 0.12345, b = 0.54321), na.action = na.omit))));
result;
}
# Stack the stacks
s <- stack(st1,st2);
# Calculate a new Raster
rNLR <- calc(s, test);
``````
• Thanks for your answer. As I said my data are temporal for one year so `st1` has 365 rasters and `st2` has 365 rasters, etc...I wonder how we change this part `x <- r[1:3]; y <- r[4:6]` ?? – usersam Mar 11 '15 at 8:24
• The model is not important but the way we do it is important.Thanks once again. – usersam Mar 11 '15 at 8:43
• If you think of r as the pixel. At that point you would have 365 layers twice in the stack. So to put the pixel value into the vector x you specify it as x <- r[1~365]. Y is then set to the next 365 layer values in the stack. – MickyT Mar 11 '15 at 8:51
• so Y would be `y <- r[366~730]`?? what If I have another variable `z`, what would be `Z <- r[?~?]` – usersam Mar 11 '15 at 8:55
• To be honest, I wasn't trying to help with the actual model, rather the way to apply it to the raster. – MickyT Mar 12 '15 at 19:53