[R-sig-Geo] How to extract standard error, r-squared for a linear model in rasterstack?

[Robert J. Hijmans]

Eddie, here is how you can do that (presented as a general worked
example for writing functions for calc) Robert

#1 create test data
library(raster)
r <- raster(nrow=10, ncol=10)
set.seed(0)
s <- stack(lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3) )))
time <- 1:nlayers(s)

#2 create an appropriate function that returns the values of interest

fun <- function(x) {
  m <- lm(x ~ time)
  s  <- summary(m)
  r2 <- s$r.squared
  p <- s$coefficients[2,4]
  cbind(r2, p)
}


#3 test the function for a single cell (vector)
test <- as.vector(s[1])
fun(test)
#            r2           p
#[1,] 0.453078 0.01643562
# OK

#4 also test with some NA values
test[3] <- NA
fun(test)
#            r2          p
#[1,] 0.4025646 0.03600484
#OK

#5 also test with only NA values
test[] <- NA
fun(test)
#Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
#  0 (non-NA) cases

#6 not good, fix function
fun <- function(x) {
if (all(is.na(x))) {
return(cbind(NA,NA))
}
m <- lm(x ~ time)
s  <- summary(m)
r2 <- s$r.squared
p <- s$coefficients[2,4]
cbind(r2, p)
}

#7 test again
test <- as.vector(s[1])
fun(test)
#            r2           p
#[1,] 0.453078 0.01643562
test[3] <- NA
fun(test)
#            r2          p
#[1,] 0.4025646 0.03600484
test[] <- NA
fun(test)
#     [,1] [,2]
#[1,]   NA   NA
# all good now

#8 test within apply (although this may not be strictly required)

t(apply(s[1:5], 1, fun))
#          [,1]         [,2]
#[1,] 0.4530780 0.0164356218
#[2,] 0.6958559 0.0007421297
#[3,] 0.6469321 0.0016099809
#[4,] 0.3184170 0.0559753606
#[5,] 0.4496749 0.0170002701
# looks good

#9 test with Raster* data
 s[1:5] <- NA
 r <- calc(s, fun)
 plot(r)

#10 all good, now try with the real data

On Sat, Nov 23, 2013 at 8:12 AM, Eddie Smith <eddieatr at gmail.com> wrote:
> Dear list,
>
> r <- raster(nrow=10, ncol=10)
> s1 <- s2<- list()
> for (i in 1:12) {
>     s1[i] <- setValues(r, rnorm(ncell(r), i, 3) )
>     s2[i] <- setValues(r, rnorm(ncell(r), i, 3) )
> }
> s1 <- stack(s1)
> s2 <- stack(s2)
>
> stest <- stack(s1, s2)
> fun <- function(x) { lm(x ~ time)$coefficients }
> x <- calc(stest, fun)
>
> >From the raster document, x contains the intercept and slope of the model.
>
> Anybody could help me in producing other details of the model e.g standard
> error, r-squared, p-value etc?
>
> Thanks in advance.
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo