-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcalculating_arima_slope_pval.R
66 lines (53 loc) · 2.33 KB
/
calculating_arima_slope_pval.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
## Running ARIMA model for the new harmonized DMSP data
## load libraries
library(raster) # for reading rasters and using parallel clusterR function
library(tictoc) # for measuring time
# Loading stack of 22 years of DMSP data cut to 45° latitude
stack_ntl <- stack("~/data/ntl/ntl_results/aurora_correction/corrected_ntl_stack.tif")
print("Start calculations")
## Calculating arima slope
fun_slope <- function(vect, time = c(1:22)) { if (sum(is.na(vect)) >= 13){ NA }
else {if (all(vect == 0)) return(0) else
ifworks = try(mod <- arima(vect, xreg=time, order = c(1, 0, 0)), silent = T)
if (class(ifworks) == "try-error") {
print(vect)
ifworks2 = try(mod <- arima(vect, xreg=time, order = c(1, 0, 0), method="ML"))
if (class(ifworks2) == "try-error") return(99)
}
return(mod$coef[3])}}
# takes approx. 2 hours on 32 cores, 32 GB RAM (16GB RAM may be enough)
tic("linear regression arima on 32 nodes")
beginCluster(n = 32)
arima_slopes = clusterR(stack_ntl, calc,
args=list(fun=fun_slope)
)
endCluster()
toc()
writeRaster(arima_slopes, "~/data/ntl/ntl_results/aurora_correction/arima_slopes.tif", overwrite=T)
## Calculating arima pvalue
fun_pval <- function(x, time = c(1:22)) {
if (((sum(is.na(x)) >= 13) | any(is.infinite(x)) | (all(x == 0)))){
NA } else {
ifworks = try(mod <- arima(x, xreg=time, order = c(1, 0, 0)), silent = T)
if (class(ifworks) == "try-error") {
ifworks2 = try(mod <- arima(x, xreg=time, order = c(1, 0, 0), method="ML"),
silent = T)
if (class(ifworks2) == "try-error") return(99)
}
pvalues = (1-pnorm(abs(mod$coef)/sqrt(diag(mod$var.coef))))*2
return(pvalues[3])
}}
tic("linear regression pval on 32 nodes")
beginCluster(n = 32)
arima_pvals = clusterR(stack_ntl, calc,
args=list(fun=fun_pval)
)
endCluster()
toc()
writeRaster(arima_pvals, "~/data/ntl/ntl_results/aurora_correction/arima_pvals.tif", overwrite=T)
# saving significant slopes
slps = raster("~/data/ntl/ntl_results/aurora_correction/arima_slopes.tif")
pvals = raster("~/data/ntl/ntl_results/aurora_correction/arima_pvals.tif")
slps[pvals>0.05]=100 # non-significant slopes are set to 100
#plot(slps)
writeRaster(slps, "~/data/ntl/ntl_results/aurora_correction/arima_slopes_significant.tif", overwrite=T)