问题描述
对于几年前回答的问题,我基本上有一个后续问题,涉及在栅格堆栈上进行线性回归。 (请参见UPDATE statement (Microsoft Access SQL)) 我进行了线性回归,并使用$ coefficients中的值计算了趋势,但现在想知道相关的p值(每个栅格像素一个)。
但是,calc抱怨无法为签名“” integer”,“ function””找到函数“ calc”的继承方法
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)
我可以理解calc不知道如何将摘要输出转换为新的栅格堆栈。因此,我尝试将lm的输出重塑为其他形式,例如使用“ broom”包,但无济于事。现在calc抱怨:is.infinite(v)中的错误:即使试图将输出强制为data.frame或as.numeric,也无法为“列表”类型实现默认方法。例如这样
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)
请帮助
解决方法
这里是一个示例,基于?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,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)