Conducting linear correlation based on stacked raster in R?

I am analyzing average annual GPP during 2000-2014 of a area. I want to conduct linear correlation analysis between annual GPP and time, and I also want to get the slope and correlation coefficient values of the linear relationship of each cell and output the results as rasters.

Do you know how to calculate the slope and R^2 values in "raster" package in R?

Can calc() or lm () will make it?

I list the data as following:

library(raster)
GPP2000=raster("GPP2000.tif")
GPP2001=raster("GPP2001.tif")
GPP2002=raster("GPP2002.tif")
GPP2003=raster("GPP2003.tif")
GPP2004=raster("GPP2004.tif")
GPP2005=raster("GPP2005.tif")
GPP2006=raster("GPP2006.tif")
GPP2007=raster("GPP2007.tif")
GPP2008=raster("GPP2008.tif")
GPP2009=raster("GPP2009.tif")
GPP2010=raster("GPP2010.tif")
GPP2011=raster("GPP2011.tif")
GPP2012=raster("GPP2012.tif")
GPP2013=raster("GPP2013.tif")
GPP2014=raster("GPP2014.tif")
GPP=stack(GPP2000,GPP2001,GPP2002,GPP20003,GPP2004,GPP2005,GPP2006,GPP2007,GPP2008,GPP20009,GPP2010,GPP2011,GPP2012,GPP2013,GPP2014)
• There's precisely an example of that in the help page of raster::calc(). See # regression of values in one brick (or stack) with 'time'. Also note that stack() can take a list of filenames as argument so that you don't need to import each layer with raster() – Loïc Dutrieux Jan 13 '17 at 9:12

Just copying the code snippet of the raster::calc() help page.

library(raster)

# Create test RasterStack
r <- raster(ncols=36, nrows=18)
r[] <- 1:ncell(r)
s <- stack(r, r*2, sqrt(r))

# regression of values in one brick (or stack) with 'time'
time <- 1:nlayers(s)
fun <- function(x) { lm(x ~ time)\$coefficients }
x2 <- calc(s, fun)

To get slope and r.squared, you'd need to adapt the function a little

fun2 <- function(x) { model <- summary(lm(x ~ time)); c(model\$coefficients[2,1], model\$r.squared) }
x2 <- calc(s, fun2)