How to extract standard error, r-squared for a linear model in rasterstack?

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

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

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20131126/04e58826/attachment.pl>