CODEWITHSUNDEEP
rlmr-rasterp-value
Emma Daniels
Emma Daniels

Reputation: 1

Extracting p-values from lineair regression on raster image through calc function (in R)

I basically have a followup question to one answered a few years ago about conducting linear regression on a raster stack. (See Linear regression on raster images - lm complains about NAs) I did the linear regression and calculated the trends with values from $coefficients but now would like to know the associated p-values (one for each raster pixel).

However, calc complains unable to find an inherited method for function ‘calc’ for signature ‘"integer", "function"’

You can reproduce this error with the following code:

library(raster)
names = c('...','...','...','...','...')
s <- stack(names)
y <- values(s)
x <- log(c(10,20,30,40,50))    
funa <- function(y) { 
        if(all(is.na(y))) {
            c(NA, NA)
        } else {
            summary(lm(y ~ x))$coefficients
        }
    }
    r <- calc(s, funa)

I can understand that calc doesn't know how the translate the output of summary to a new raster stack. So I've tried reshaping the output of lm to other forms, for example with the "broom" package, to no avail however. Now calc complains that: Error in is.infinite(v) : default method not implemented for type 'list' even when I try to force the output to a data.frame or as.numeric. Like this for example

library(broom)
funlm <- function(y) {
    if(all(is.na(y))) {
        c(NA, NA)
    } else {
        as.data.frame(glance(lm(y ~ x)))
    }
}
r <- calc(s, funlm)

Please help

Upvotes: 0

Views: 272

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47146

Here is an example, based on the example in ?calc

# create data
library(raster)
r <- raster(nrow=10, ncol=10)
s1 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))
s2 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))
s1 <- stack(s1)
s2 <- stack(s2)

# regression of values in one brick (or stack) with another, extract p-values
s <- stack(s1, s2)
fun <- function(x) {
    if (all(is.na(x))) {
       return(c(NA, NA))
    } 
    m <- lm(x[1:12] ~ x[13:24])
    summary(m)$coefficients[,4]
}
x1 <- calc(s, fun)

Upvotes: 0

Related Questions