Model-average predictions for a habitat.stack (raster)

[Patrick Johnson]

Hey everybody - what is the best way to make maps using model-averaged predictions in unmarked? I have had success making model-averaged predictions using a fitList on a dataframe (as spelled out in the reference manual), but when I tried using similar methods to apply to a raster layer (a RasterStack) I get the error: 

"Error in x[, "Predicted"] : object of type 'S4' is not subsettable" 

My models with the greatest AIC support include 4-6 environmental variables. 

Any advice would be appreciated, thanks! 

[Richard Chandler]

Hi Pat,

This hasn't been implemented yet. I'll put it on the to-do list. In the meantime, you can convert your raster data into a data.frame (rows are pixels, columns are covariate) and then supply that data.frame to predict(). 

Richard

[Patrick Johnson]

Thanks Richard - good idea on turning the raster into a data.frame, I hadn't thought of that. Thanks! 

[Chris Hallam]

Hi Patrick,

I found your post here and was having a similar issue.  I can predict to the dataframe when I convert the stack.  I was wondering if there is a way of turing the predicitons back into a raster for plotting purposes?

I realize its a long time ago but though Id ask...

Cheers

Chris

[Patrick Johnson]

Chris, 

It has been a while for me, so maybe someone else can jump in and help if this doesn't. I think you can use the writeRaster function in the Raster Package (http://www.inside-r.org/packages/cran/raster/docs/writeRaster) to convert your data back to a raster. Then you can plot the data in Arc or some other GIS. 

Hopefully this helps and if not I hope you find out what you need. Good luck with your work. 

Best,

Pat  

--
You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/q3FO84fbyGI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Chris Sutherland]

rasterFromXYZ () in the raster package also does this type of thing. You simply provide a matrix with 3 columns - the xy coordinates and then the predictions at each location.

You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.

[Chris Hallam]

Hi Thanks,

I managed to plot the model averaged results by running all the models seperately and making a stack and then taking the mean of the predictions...I hope that would be ok?

I converted the results to a data frame and can extract estimates of N etc (Im running occuRN)

Thanks for your help all...

Chris

[Chris Sutherland]

If i understand what you have done correctly, i think this should be a _weighted_ mean (I.e. multiplying the model predictions by the model weights).

[Norberello A]

Hi Richard,

I am aware this post is 5yr old and perhaps is dead, but has this been implemented already? I've got exactly the same problem today and if possible I don't want to go over the data.frame solution.

Thanks for the answer in advance,

N

[Richard Chandler]

Hi Norberello,

I think should work. Have you tried it? If it doesn't work, you can cover your raster stack to a data.frame like this:

mydat <- as.data.frame(myrast, xy=TRUE)

Richard

--

You received this message because you are subscribed to the Google Groups "unmarked" group.

To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Richard Chandler

Assistant Professor

Warnell School of Forestry and Natural Resources

University of Georgia

706-542-5815

chandlerlab.uga.edu

[Norberello A]

Thanks for the quick reply, Richard,

I see, then it should a problem of any of the rasters not having the same # of pixels or something of that sort. After all, it was working fine, it only gave that error when only about 700 predicition lines out of 20,000  were left to finish. Still, when I check if the rasters used for stack are equally OK by

compareRaster(sas.p, Lar.p, For.in.p,HQ.p,elev.p,

              stream.p,slope.p)

everything seems correct, and I had not problems after when putting them together in the stack.

Well, I ended up with the dataframe solution, perhaps not in the most elegant way as I didn't manage to transform the raster stack into a dataframe in only one step: 

habitat <- stack(sas.p, Lar.p, For.in.p,HQ.p,elev.p,

                 stream.p,slope.p)

names(habitat) <- c("sas.p", "Lar", "For.in","HQ","Elev",

                    "Stream","Slope")

habitat.t<-data.frame(rasterToPoints(habitat))

#this didn't work, too many NAs, I doesn't make sense

habitat.t2<-lapply(habitat,as.data.frame) 

#same problem

habitat.t3 <- as.data.frame(habitat, xy=TRUE) 

#only get x and y, but the covariate columns are all NAs :(

So, I converted every raster into a dataframe first which I merged afterwards to be used as the predictor table:

sas.p.t<-data.frame(rasterToPoints(sas.p))

lar.p.t<-data.frame(rasterToPoints(Lar.p))

for.in.t<-data.frame(rasterToPoints(For.in.p))

HQ.p.t<-data.frame(rasterToPoints(HQ.p))

elev.p.t<-data.frame(rasterToPoints(elev.p))

stream.p.t<-data.frame(rasterToPoints(stream.p))

slope.p.t<-data.frame(rasterToPoints(slope.p))

#and now to put them together, 

##I know, I know, this is messy, but I couldn't figure a better wa

mytempdata = merge(sas.p.t, lar.p.t,by=c("x","y"))

mytempdata = merge(mytempdata,for.in.t,by=c("x","y"))

mytempdata = merge(mytempdata,HQ.p.t,by=c("x","y"))

mytempdata = merge(mytempdata,elev.p.t,by=c("x","y"))

mytempdata = merge(mytempdata,stream.p.t,by=c("x","y"))

mytempdata = merge(mytempdata,slope.p.t,by=c("x","y"))

names(mytempdata) <- c("x","y","sas.p", "Lar", "For.in",

                       "HQ","Elev",

                    "Stream","Slope")

write.csv(mytempdata,file="dataframe.csv")

I still needed to set the dataframe properly in excel before using it for the predicition by saving it to a file first.

write.csv(mytempdata,file="dataframe.csv")

these are the averaged models by the way:

fmList <- fitList(fm95,fm27,fm55,fm62,fm45,fm8)

and so I could run the prediction:

E3 <- predict(fmList, type="state", newdata=mytempdata)

head(E3)

#good!

#save it as csv so you can work in excel to paste X and Y

write.csv(E3,file="E3.csv")

This didn't end here as the prediction doesn't come with X and Y (is there an automatic way to do so? Can you add into the output predicted, SD, Max, Min columns the two extra columns with the corresponding X and Y?).

So, anywasys, I used the X and Y back from the mytempdata dataframe. I wasn't clever enough to do it in R and I just combined the two files in excel, save as csv, and entered in QGIS for the map (attached).

It worked very well, but probably you think I am crazy for doing so much back an forth. Any advice to improve these steps very welcome,

Norber

write.csv(E3,file="E3.csv")

On Wednesday, 7 June 2017 22:20:56 UTC+7, Richard Chandler wrote:

Hi Norberello,

I think should work. Have you tried it? If it doesn't work, you can cover your raster stack to a data.frame like this:

mydat <- as.data.frame(myrast, xy=TRUE)

Richard

On Wed, Jun 7, 2017 at 2:45 AM, Norberello A <norbe...@gmail.com> wrote:

Hi Richard,

I am aware this post is 5yr old and perhaps is dead, but has this been implemented already? I've got exactly the same problem today and if possible I don't want to go over the data.frame solution.

Thanks for the answer in advance,

N

On Monday, 29 October 2012 20:04:52 UTC+7, Richard Chandler wrote:

Hi Pat,

This hasn't been implemented yet. I'll put it on the to-do list. In the meantime, you can convert your raster data into a data.frame (rows are pixels, columns are covariate) and then supply that data.frame to predict(). 

Richard

 

On Sat, Oct 27, 2012 at 2:46 PM, Patrick Johnson <patrickly...@gmail.com> wrote:

Hey everybody - what is the best way to make maps using model-averaged predictions in unmarked? I have had success making model-averaged predictions using a fitList on a dataframe (as spelled out in the reference manual), but when I tried using similar methods to apply to a raster layer (a RasterStack) I get the error: 

"Error in x[, "Predicted"] : object of type 'S4' is not subsettable" 

My models with the greatest AIC support include 4-6 environmental variables. 

Any advice would be appreciated, thanks! 

--

You received this message because you are subscribed to the Google Groups "unmarked" group.

To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Richard Chandler]

What error message did you get when it failed with just 700 rows left? It could be a memory issue. Also, you can get the x and y coordinates using: 

dat <- as.data.frame(rasters, xy=TRUE)

To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Norberello A]

Thanks, Richard! I will try it in a different computer, let you know if all runs properly,

cheers,
Norber

You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/q3FO84fbyGI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+u...@googlegroups.com.

[Norberello A]

Ups, sorry, forgot to answer, it was the same error message as above:

"Error in x[, "Predicted"] : object of type 'S4' is not subsettable”

N

On 9 Jun 2017, at 07:43, Richard Chandler <rcha...@warnell.uga.edu> wrote:

You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/q3FO84fbyGI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+u...@googlegroups.com.

[Marconi Campos Cerqueira]

Hello group,

I am aware that this post is old and perhaps is dead, but I`m getting the same error described by Norberello A

I saw an updated version (2019) of the manual to create predictions using rasters, but  this is what I get:

doing row 15000 of 16324 

  doing row 16000 of 16324 

Error in x[, "Predicted"] : object of type 'S4' is not subsettable

Thanks for comments in advance,

Best,

Marconi

To unsubscribe from this group and stop receiving emails from it, send an email to unma...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.



--
Richard Chandler
Assistant Professor
Warnell School of Forestry and Natural Resources
University of Georgia
706-542-5815
chandlerlab.uga.edu

--
You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/q3FO84fbyGI/unsubscribe.

To unsubscribe from this group and all its topics, send an email to unma...@googlegroups.com.

[Ken Kellner]

Hi Marconi,

I think you are seeing that error because you are trying to generated a model-averaged prediction using a raster as the newdata. This is still currently not supported in unmarked. If this is the case you can try converting your raster to a data frame as Richard suggests. 

Alternatively I wrote up a quick function that should allow for predicting from fitlists with a raster, while still getting a raster output. See example below based on Richard's spatial vignette. It seems to work for me - I'll write a patch to get it into the next version of unmarked so this is handled automatically.

Ken

Example:

library(unmarked)

library(raster)

data(crossbill)

data(Switzerland)

#Fit two models 

umf <- unmarkedFrameOccu(

y=as.matrix(crossbill[,c("det991", "det992", "det993")]),

siteCovs=crossbill[,c("ele", "forest")],

obsCovs=list(date=crossbill[,c("date991", "date992", "date993")]))

sc <- scale(siteCovs(umf))

siteCovs(umf) <- sc

fm.occu <- occu(~date ~ele + I(ele^2) + forest, umf)

fm.occu2 <- occu(~1~1, umf)

#Build newdata raster stack

elevation <- rasterFromXYZ(Switzerland[,c("x","y","elevation")],

    crs="+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")

forest <- rasterFromXYZ(Switzerland[,c("x","y","forest")],

    crs="+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")

ele.s <- (elevation-1189)/640

forest.s <- (forest-34.7)/27.7

ef <- stack(ele.s, forest.s)

names(ef) <- c("ele", "forest")

#Make fitlist

fl = fitList(fm.occu,fm.occu2)

#New function for model-averaged prediction from a fitList using a raster input

pred_raster_fitlist <- function(object, type, newdata=NULL,

    backTransform = TRUE, appendData = FALSE, level=0.95){

  

  #----Same as existing function

  fl <- object@fits

  ese <- lapply(fl, predict, type = type, newdata = newdata,

                backTransform = backTransform, level=level)

  #---------------------------------

  if(class(newdata) == "RasterStack"){

    if(!require(raster)) stop("raster package is required")

    ese <- lapply(ese, as.matrix)

  }

  

  #----Same as existing function

  E <- sapply(ese, function(x) x[,"Predicted"])

  SE <- sapply(ese, function(x) x[,"SE"])

  lower <- sapply(ese, function(x) x[,"lower"])

  upper <- sapply(ese, function(x) x[,"upper"])

  ic <- sapply(fl, slot, "AIC")

  deltaic <- ic - min(ic)

  wts <- exp(-deltaic / 2)

  wts <- wts / sum(wts)

  parav <- as.numeric(E %*% wts)

  seav <- as.numeric(sqrt(SE^2 + (E - parav)^2) %*% wts)

  out <- data.frame(Predicted = parav, SE = seav)

  out$lower <- as.numeric(lower %*% wts)

  out$upper <- as.numeric(upper %*% wts)

  #------------------------------

  if(class(newdata) == "RasterStack"){

    E.mat <- matrix(out[,1], dim(newdata)[1], dim(newdata)[2], byrow=TRUE)

    E.raster <- raster::raster(E.mat)

    raster::extent(E.raster) <- raster::extent(newdata)

    out.rasters <- list(E.raster)

    for(i in 2:ncol(out)) {

      i.mat <- matrix(out[,i], dim(newdata)[1], dim(newdata)[2], byrow=TRUE)

      i.raster <- raster::raster(i.mat)

      raster::extent(i.raster) <- raster::extent(newdata)

      out.rasters[[i]] <- i.raster

    }

    out.stack <- stack(out.rasters)

    names(out.stack) <- colnames(out)

    raster::crs(out.stack) <- raster::crs(newdata)

    return(out.stack)

  }

}

#Generate model-averaged prediction

raster_test <- pred_raster_fitlist(fl, type='state', newdata=ef)

plot(raster_test)

[Marconi Campos Cerqueira]

Hello Ken,

Thank you very much for your quick reply.

I was indeed trying to generate a model-averaged prediction using a raster stack as the new data.

The new function worked perfectly. Thank you very much!

Best

Marconi

[Ken Kellner]

Great! Here's the patch for future reference. 

https://github.com/rbchan/unmarked/pull/140