wish: multivariate stats command

[wegmann]

a command for linear regression of x~y1+y2+... raster would be handy and optional output of residuals (maybe optional with piecewise regressions, robust linear regressions etc.). To ease this code snippet:

   library(raster)
   library(RStoolbox)
   
   lsat88 <- brick("p224r63_1988.grd")
   lsat11 <- brick("p224r63_2011.grd")
   
   lsat88.ndvi <- spectralIndices(lsat88, red="B3_sre",nir="B4_sre",indices="ndvi")
   lsat11.ndvi <- spectralIndices(lsat11, red="B3_sre",nir="B4_sre",indices="ndvi")

   lm.lsat <- lm(lsat88.ndvi[]~lsat11.ndvi[], na.action = NULL) # can also be done as multivariate regression: x~y1+y2+y3 with dependent and explanatory variables
   summary(lm.xy)

   Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
   yc.ndvi[]   0.248039   0.001339   185.3   <2e-16
   ...
   Multiple R-squared:  0.1702,	Adjusted R-squared:  0.1702 
   ...
   
   resid.lsat <- raster(lsat88)	
   resid.lsat[] <- lm.lsat$residuals
   plot(resid.lsat,stretch="lin")