2013-11-28 47 views
2

我正在使用R來執行我的迴歸。我成功地將r平方剩餘標準錯誤提取到了rasterbrick中。然後我需要運行另一個代碼來獲得p值(F stat)。我怎麼能把fun1和fun2結合起來,這樣我就可以製作一個包含這些信息的rasterbrick?如何將回歸摘要(例如p值和coeff)輸出到rasterbrick中?

這裏是我的代碼:

library(raster) 

#1 create test data 
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) 
s[1:5] <- NA 

#2 Run function1 to obtain r-squared and residual standard error 
fun1 <- function(x) { 
if (all(is.na(x))) { 
return(cbind(NA,NA)) 
} 
m = lm(x~time) 
s <- summary(m) 
r2 <- s$r.squared 
resid.s.e <- s$sigma 
cbind(r2, resid.s.e) 
} 

#3 Run function to obtaion p-value(from F stat) 
fun2 <- function(x) { 
if (all(is.na(x))) { 
return(cbind(NA,NA)) 
} 
m = lm(x~time) 
s <- summary(m) 
r2 <- s$r.squared 
pf<- pf(s$fstatistic[1], s$fstatistic[2], s$fstatistic[3],lower.tail = FALSE) 
cbind(r2, pf) 
} 

#Apply both functions with rasterstack and plot 
r <- calc(s, fun) 
plot(r) 

r2 <- calc(s, fun2) 
plot(r2) 

在此先感謝。

+0

你的文字說p值,你的標題說係數。這是什麼? –

+1

謝謝@羅曼·盧斯特里克。編輯我的標題和帖子。我很匆忙,混了起來。道歉。 – Eddie

+1

爲什麼不讓1個函數返回'resid.s.e'和'pf',然後選擇變量傳遞給'calc'?重新計算相同的'lm'似乎是多餘的。 – zx8754

回答

2

我想我得到了答案。

在cbind()中添加更多的列將允許我在輸出rasterstack中添加更多層。

library(raster) 

#1 create test data 
    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) 
    s[1:5] <- NA 

#2 Run function1 to obtain r-squared, residual standard error and p-value(F stat) 
    fun <- function(x) { 
    if (all(is.na(x))) { 
    return(cbind(NA,NA,NA)) 
    } 
    m = lm(x~time) 
    s <- summary(m) 
    r2 <- s$r.squared 
    resid.s.e <- s$sigma 
    pf<- pf(s$fstatistic[1], s$fstatistic[2], s$fstatistic[3],lower.tail = FALSE) 
    cbind(r2, resid.s.e, pf) 
    } 

    r <- calc(s, fun) 
    r 
    class  : RasterBrick 
    dimensions : 10, 10, 100, 3 (nrow, ncol, ncell, nlayers) 
    resolution : 36, 18 (x, y) 
    extent  : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) 
    coord. ref. : +proj=longlat +datum=WGS84 
    data source : in memory 
    names  :  layer.1,  layer.2,  layer.3 
    min values : 1.300285e-01, 1.457297e+00, 5.199987e-07 
    max values : 0.9271788, 5.0219805, 0.2495312 
相關問題