I am performing a regression analysis with three raster stacks in R. This worked so far. Now I want to check the conditions of the model and I am not sure how exactly that should work in this case. Since the model is calculated individually for each pixel, I cannot extract the residues for the entire investigation area and then check them. Therefore I currently have to extract the data of the variables for a few points/coordinates on a trial basis and then create a linear model (lm) on which apply the tests (Shapiro-test, Breusch-Pagan-test, and VIF-test).

Does anyone know a more elegant solution?

NAO_stack <- stack("NAO_stack_fertig.tif")
AOD_stack <- stack("AOD_stack_fertig.tif")
CFC_stack <- stack("CFC_stack_fertig.tif")
NAO_AOD_stack <- stack(NAO_stack, AOD_stack)
NAO_CFC_stack <- stack(NAO_stack, CFC_stack)

NAO_AOD_CFC <- stack(NAO_stack, AOD_stack, CFC_stack)
NAO_CFC_AOD <- stack(NAO_stack, CFC_stack, AOD_stack)

fun1=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]);summary(m)$coefficients[,1] }}
Steigung_CFC <- calc(NAO_CFC_AOD, fun1)
names(Steigung_CFC) <- c("Intercept", "NAO_Steigung", "AOD_Steigung")

fun2=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]);summary(m)$coefficients[,4] }}
p_value_CFC <- calc(NAO_CFC_AOD, fun2)
names(p_value_CFC) <- c("Intercept", "NAO_P_Value", "AOD_P_Value")

### pixel-by-pixel multiple linear regression with NAO and AOD / CFC as dependent variables

fun3=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]);summary(m)$adj.r.squared }}
r_squared_CFC <- calc(NAO_CFC_AOD, fun3)
plot(r_squared_CFC, main = "R²_CFC_multiple")

r_squared_AOD <- calc(NAO_AOD_CFC, fun3)
plot(r_squared_AOD, main = "R²_AOD_multiple")

### pixel-wise regression only with NAO as a dependent variable

fun4=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55]);summary(m)$adj.r.squared }}
r_squared_CFC_NAO <- calc(NAO_CFC_stack, fun4)
plot(r_squared_CFC_NAO, main = "R²_CFC")

r_squared_AOD_NAO <- calc(NAO_AOD_stack, fun4)
plot(r_squared_AOD_NAO, main = "R²_AOD")

r_squared_AOD_CFC <- calc(stack(CFC_stack, AOD_stack), fun4)
plot(r_squared_AOD_CFC, main = "R²_AOD_CFC")

### Extract data for a point to check the assumptions on the resulting LM model
XCoordinate = 8
YCoordinate = 61
points(XCoordinate, YCoordinate, pch = 15, cex = 0.5)
xy <- cbind(XCoordinate,YCoordinate)
sp <- SpatialPoints(xy)
data_CFC <- extract(CFC_stack, sp)  
data_CFC <- as.vector(data_CFC)
data_AOD <- extract(AOD_stack, sp) 
data_AOD <- as.vector(data_AOD)
data_NAO <- extract(NAO_stack, sp)
data_NAO <- as.vector(data_NAO)

CFC_lm <- lm(data_CFC ~ data_AOD + data_NAO)
CFC_resid <- CFC_lm$residuals

AOD_lm <- lm(data_AOD ~ data_CFC + data_NAO)
AOD_resid <- AOD_lm$residuals

  • Well, you do not have a regression solution for the entire raster so, none of the aforementioned statistics are relevant. Each pixel is an independent regression so things like Shapiro, Breusch-Pagan and VIF are only relevant at the pixel-level as well. I do have some issues with this type of raster regression approach. Why not just take a point sample of the rasters, specify a regression, evaluate competing model, evaluate fit and the predict the final model to the rasters? This is more supported and lets you know if your model is valid or not. Commented Jun 2, 2020 at 23:22

1 Answer 1


ive solved my problem with the following code:

## Review of the assumptions of the multiple linear regression models

### Shapiro-Test
fun5=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]); shapiro.test(m$residuals)$p.value}}
nv_CFC_lm <- calc(NAO_CFC_AOD, fun5)
mask_nv_CFC_lm <- nv_CFC_lm
mask_nv_CFC_lm[mask_nv_CFC_lm > 0.05] <- NA

nv_AOD_lm <- calc(NAO_AOD_CFC, fun5)
mask_nv_AOD_lm <- nv_AOD_lm
mask_nv_AOD_lm[mask_nv_AOD_lm > 0.05] <- NA

### BP-Test 
fun5_1=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]); bptest(m)$p.value}}
bp_CFC_lm <- calc(NAO_CFC_AOD, fun5_1)
mask_bp_CFC_lm <- bp_CFC_lm
mask_bp_CFC_lm[mask_bp_CFC_lm > 0.05] <- NA

bp_AOD_lm <- calc(NAO_AOD_CFC, fun5_1)
mask_bp_AOD_lm <- bp_AOD_lm
mask_bp_AOD_lm[mask_bp_AOD_lm > 0.05] <- NA

### VIF-Test
fun5_2=function(x) { if (is.na(x[1])){ NA } else { m <- lm(x[56:110] ~ x[1:55] + x[111:165]); vif(m)}}
vif_CFC_lm <- calc(NAO_CFC_AOD, fun5_2)

vif_AOD_lm <- calc(NAO_AOD_CFC, fun5_2)

The output is an raster with the p-Value (Shapiro and BP-Test) or the VIF-Value for my study area. So i can detect areas, where the conditions of the model are not met.

@JeffreyAdams so you mean, that you would for every point of the rasterstacks the values and would build a linear model with those values for every coordinate-pair?

  • No, I mean that you should take some sort of sample and run a global model. Functionally, you are dealing with a population and it really should be a sample. The other thing, I cannot imagine what you are representing but it looks like, for each pixel, you are drawing n=55 values for y and then 55 values for each raster stack. How can this possibly be independent and not severely pseudoreplicated? If you can figure out how to specify a global model you can at lease evaluate the validly of said model. Commented Jun 3, 2020 at 13:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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