T-Test for population raster in R

Question

I am trying to perform a "t-test" for 4 population rasters for the entire tropical strip for my graduate thesis. Sample years are 2000-2005-2010-2015. I am using this script:

setwd("C:/Users/MICITT/Desktop/TESIS_MSC_MARCOMM/Poblacion")
require(raster)
require(rgdal)
pob00 <- raster("PRUEBA_rem_pob00_SC.tif")
pob05 <- raster("PRUEBA_rem_pob05_SC.tif")
pob10 <- raster("PRUEBA_rem_pob10_SC.tif")
pob15 <- raster("PRUEBA_rem_pob15_SC.tif")
pob_ttest <- t.test(pob00, pob05, pob10, pob15, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95)
writeRaster(pob_ttest, filename = "C:/Users/MICITT/Desktop/TESIS_MSC_MARCOMM/Poblacion/ttestprueba", format="GTiff", overwrite=TRUE"

However, when I try to save the raster, it displays the following error:

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘writeRaster’ for signature ‘"htest", "character"’

I need to get a final raster that gives me the value of "t", "df" and "p-value" per pixel. Anyone know how to solve this?

Answer 1

you should conceptually review what you really want to test across the 4 rasters. The only way I see so that the output of the test is a raster object itself is that you want to test whether the unit in each pixel is statistically different from zero, by sampling the values of each pixel across the 4 rasters. If that is the case, you can try the following (very clumsy but illustrates the idea):

require(raster)
require(rgdal)
require(data.table)
rs <- raster::stack(pob00, pob05, pob10, pob15)
df.rs <- data.table(values(rs))
# adding the variables you want to calculate:
df.rs$df <- df.rs$t <- df.rs$pvalue <- as.numeric(NA)
# now you do the t-test in each pixel (can take very long if rasters have numerous pixels):
for (i in 1:nrow(df.rs)) {
  test <- t.test(as.numeric(df.rs[i,-c('df','t','pvalue')]))
  df.rs[i,'df'] <- as.numeric(test$parameter)
  df.rs[i,'t'] <- as.numeric(test$statistic)
  df.rs[i,'pvalue'] <- as.numeric(test$p.value)
}
# creating the outcomes as a raster stack (each layer will be a statistic):
r.t <- raster::stack(pob00,pob00,pob00)
values(r.t) <- as.matrix(df.rs[,c('pvalue', 't', 'df')])
# this should work for you; to plot them:
plot(r.t[[1]]) # p-value
plot(r.t[[2]]) # t-stat
plot(r.t[[3]]) # df