

Hi,
I have a GAM model that contains factors in addition to smooth terms etc. I
would like to evaluate this at all points on a raster, but I can't seem to
get the predict.raster() function to pick up the values of my factors. I
have created a minimum (not) working example here, based on the example for
the predict() function.
Any suggestions in this direction would be appreciated.
Mark
# create a RasterStack or RasterBrick with with a set of predictor layers
logo < brick(system.file("external/rlogo.grd", package="raster"))
# known presence and absence points
p < matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85,
66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46, 38,
31,
22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)
a < matrix(c(22, 33, 64, 85, 92, 94, 59, 27, 30, 64, 60, 33, 31, 9,
99, 67, 15, 5, 4, 30, 8, 37, 42, 27, 19, 69, 60, 73, 3, 5, 21,
37, 52, 70, 74, 9, 13, 4, 17, 47), ncol=2)
# extract values for points
xy < rbind(cbind(1, p), cbind(0, a))
v < data.frame(cbind(xy[,1], extract(logo, xy[,2:3])))
colnames(v)[1] < 'pa'
#Convert blue layer to a factor
v$blue < factor(v$blue < 200)
#build a model, here an example with glm
model < glm(formula=pa~., data=v)
#Setup prediction objects
pred.b < logo[[1:2]]
const.df < data.frame(blue=TRUE)
#predict to a raster
r1 < predict(pred.b, model, progress='text',
const=const.df)
#Returns an error as follows:
# Error in `[.data.frame`(blockvals, , f[j]) : undefined columns selected
[[alternative HTML version deleted]]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


Mark,
I had not considered the case where the constant is a factor. I need
to fix that. Here is a script that accomplishes what you want, I think
# known presence and absence points
p < matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85,
66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46, 38,31,
22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)
a < matrix(c(22, 33, 64, 85, 92, 94, 59, 27, 30, 64, 60, 33, 31, 9,
99, 67, 15, 5, 4, 30, 8, 37, 42, 27, 19, 69, 60, 73, 3, 5, 21,
37, 52, 70, 74, 9, 13, 4, 17, 47), ncol=2)
# extract values for points
xy < rbind(cbind(1, p), cbind(0, a))
v < data.frame(cbind(xy[,1], extract(logo, xy[,2:3])))
colnames(v)[1] < 'pa'
#Convert blue layer to a factor
## RH I would not use TRUE/FALSE but an integer representation instead
v$blue < factor((v$blue < 200) + 0)
#build a model, here an example with glm
model < glm(formula=pa~., data=v)
#Setup prediction objects
pred.b < logo
pred.b[[3]] < (pred.b[[3]] < 200) + 0
# or all TRUE? pred.b[[3]] < 1
names(pred.b) < names(logo)
#predict to a raster
r1 < predict(pred.b, model, progress='text')
You say you want to use a gam, but the example is a glm. With gam
there might be additional issues. Please let me know if that is the
case.
Robert
On Mon, Apr 20, 2015 at 3:47 AM, Mark Payne < [hidden email]> wrote:
> Hi,
>
> I have a GAM model that contains factors in addition to smooth terms etc. I
> would like to evaluate this at all points on a raster, but I can't seem to
> get the predict.raster() function to pick up the values of my factors. I
> have created a minimum (not) working example here, based on the example for
> the predict() function.
>
> Any suggestions in this direction would be appreciated.
>
> Mark
>
> # create a RasterStack or RasterBrick with with a set of predictor layers
> logo < brick(system.file("external/rlogo.grd", package="raster"))
>
> # known presence and absence points
> p < matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85,
> 66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46, 38,
> 31,
> 22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)
>
> a < matrix(c(22, 33, 64, 85, 92, 94, 59, 27, 30, 64, 60, 33, 31, 9,
> 99, 67, 15, 5, 4, 30, 8, 37, 42, 27, 19, 69, 60, 73, 3, 5, 21,
> 37, 52, 70, 74, 9, 13, 4, 17, 47), ncol=2)
>
> # extract values for points
> xy < rbind(cbind(1, p), cbind(0, a))
> v < data.frame(cbind(xy[,1], extract(logo, xy[,2:3])))
> colnames(v)[1] < 'pa'
>
> #Convert blue layer to a factor
> v$blue < factor(v$blue < 200)
>
> #build a model, here an example with glm
> model < glm(formula=pa~., data=v)
>
> #Setup prediction objects
> pred.b < logo[[1:2]]
> const.df < data.frame(blue=TRUE)
>
> #predict to a raster
> r1 < predict(pred.b, model, progress='text',
> const=const.df)
>
>
> #Returns an error as follows:
> # Error in `[.data.frame`(blockvals, , f[j]) : undefined columns selected
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo

