N-mixture model maps for abundace: interpolate covariates as rasters?

[Nicole Angeli]

Hello all;

I am modeling lizard abundance in 1260 sq m patches across a small island. Each of 41 sites was sampled 6 times (double independent observers). I am using each observer's survey as a site visit, and we collected 12 site (environmental) and 1 survey (observer) covariate. After using the dredge function on models for abundance (~1 ~..covariates..) or detection (~.covariates.. ~1), I used the top abundance and detection covariates in a model. The model doesn't improved when I add or delete covariates on either side. I used parboot with 100 simulations and GOF was fine. The model is here:

> mAbunD

Call:
pcount(formula = ~Distance.from.Enclosure + Cover.Percent + Habitat +
    AvgT.OTM + MaxT.Hobo + MaxT.OTM ~ Distance.from.Enclosure +
    LLDepth + Woody.stems + AvgT.Hobo + AvgT.OTM, data = buckRM)

Abundance:
                                        Estimate       SE      z     P(>|z|)
(Intercept)                         4.37476 2.549431   1.72   8.62e-02
Distance.from.Enclosure   -0.00473 0.000373 -12.67   8.44e-37
LLDepth                          -0.19061 0.041985  -4.54    5.63e-06
Woody.stems                 -0.05090 0.016867  -3.02     2.55e-03
AvgT.Hobo                     -0.24658 0.102684  -2.40    1.63e-02
AvgT.OTM                       0.26193 0.087816   2.98    2.86e-03

Detection:
                                           Estimate       SE     z     P(>|z|)
(Intercept)                          3.969698 2.172067  1.83   0.067608
Distance.from.Enclosure   -0.000982 0.000283 -3.48    0.000511
Cover.Percent                      -0.009896 0.002592 -3.82 0.000135
Habitat2                            0.250099 0.113097  2.21    0.027011
Habitat3                          -0.279203 0.099127 -2.82    0.004853
AvgT.OTM                         -0.270873 0.089419 -3.03   0.002452
MaxT.Hobo                        0.021481 0.005779  3.72   0.000201
MaxT.OTM                         0.037335 0.010725  3.48  0.000499

AIC: 910.0409

Then, I found site abundance using:

>abund<-ranef(mAbunD, K=500)
>EBUP<-bup(abund,stat="mode")
> sum(bup(abund))  ##get abundance estimate
>CI<-confint(abund, level=0.9)
> colSums(CI)


I am plotting curves of the individual covariates using this sort of code:
>newData2 <- data.frame(AvgT.OTM=seq(25,50,by=0.5), Cover.Percent=0, Distance.from.Enclosure=0,Habitat=factor(2,levels=c(1,2,3)), MaxT.Hobo=0,MaxT.OTM=0)
>E.p <- predict(mAbunD, type="state", newdata=newData2, appendData=TRUE)


The next step is estimating the total abundance across the island.

Based on other posts on this listserv (e.g.) I am wondering if it would be appropriate to interpolate the covariates directly on the island (in GIS or R) and then stack those interpolated surfaces as rasters to map abundance with 'predict' (see vignette)?

I appreciate any insights into my process here...

Please let me know if you need any further information to help.

Best,
Nicole

[Jeffrey Royle]

hi Nicole,

 That may be reasonable for smoothly varying covariates but possibly not so much for covariates that are very patchy (i.e., without much spatial corrleation). That said, it might be  your only option if you don't have the GIS layers for the covariate information.

 We have a couple examples of this type of prediction problem on the unmarked Google Site, under the PWRC workshop materials.  There is an example using hierarchical distance sampling (the Island Scrub Jay) and another example using a multinomial type of model applied to mapping bird abundance from the Swiss Breeding Bird Survey. Check those examples out if you run into any problems.

 

regards

andy

 

--
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/groups/opt_out.

[Nicole Angeli]

Hello;

Thanks so much! I was able to rasterize data by calibrating it against derived Landsat temperatures, SSURGO data, etc., and then using those rasters in the pcount models with similar results to what I found from the collected data.

However, now I try to use predict after successfully stacking the rasters:

#mAbunDA<- pcount(formula = ~Sand.Content + Habitat + Elevation + Water.Content + NDVI + Temperature ~ Sand.Content + Habitat + Elevation + Water.Content + NDVI + Temperature, data = buckRM)

Abundance:
              Estimate       SE      z  P(>|z|)

(Intercept)   -10.5497  2.32254 -4.542 5.56e-06
Sand.Content   -2.4227  0.33402 -7.253 4.07e-13
Habitat3       -0.1110  0.17917 -0.619 5.36e-01
Habitat4       -1.3209  0.21153 -6.244 4.25e-10
Habitat5       -6.4687 14.86247 -0.435 6.63e-01
Habitat6       -0.2299  0.42085 -0.546 5.85e-01
Elevation      -0.0248  0.00383 -6.475 9.47e-11
Water.Content   2.4137  0.34879  6.920 4.51e-12
NDVI            2.0060  1.43372  1.399 1.62e-01
Temperature     0.6085  0.10932  5.567 2.60e-08

Detection:
              Estimate      SE       z P(>|z|)

(Intercept)   -2.05564 1.75385 -1.1721  0.2412
Sand.Content  -0.01309 0.28245 -0.0464  0.9630
Habitat3      -0.12639 0.15191 -0.8320  0.4054
Habitat4       0.12297 0.17446  0.7049  0.4809
Habitat5      -0.47688 0.33287 -1.4326  0.1520
Habitat6       0.64399 0.43503  1.4803  0.1388
Elevation     -0.00200 0.00246 -0.8097  0.4181
Water.Content -0.04317 0.28870 -0.1495  0.8811
NDVI           2.47593 1.13962  2.1726  0.0298
Temperature    0.00629 0.08514  0.0739  0.9411

AIC: 991.7011


pred.data<-stack(step1, step2, step3, sand.crop,step4, step5)  ##use rasters 'stepX' resampled to sand.crop to force extent
names(pred.data) <- c("Water.Content", "Elevation", "Habitat", "Sand.Content","NDVI", "Temperature")

E <- predict(mAbunDA, type="state", newdata=pred.data)
Error: ncol(coefficients) == length(obj@estimates) is not TRUE

When I look at the data,

#str(pred.data)
Formal class 'RasterStack' [package "raster"] with 11 slots

#str(mAbunDA)
Formal class 'unmarkedFitPCount' [package "unmarked"] with 14 slots

I assume that because 'Habitat' is a factor with 6 levels, I am running into problems, although I don't necessarily know if that's true. I can 'Extract Multivalue Pts' in ArcGis and successfully run the predict function for each of the six levels using a dataframe, but it results in running the models using 'predict'  for each variable for each habitat type which isn't too useful. Maybe that's the way to go, though.

I could send you reproducible data plus rasters, sorry I haven't done that here. Any insights are appreciated!

Nicole

[Richard Chandler]

Hi Nicole,

I don't think the predict function can handle raster stacks that include factors. However, I think you could convert your factors to dummy variables and use the dummy variables for model fitting and prediction. 

Richard

--

Richard Chandler

Assistant Professor

Warnell School of Forestry and Natural Resources

University of Georgia

706-542-5815

chandlerlab.wordpress.com

[Chris Sutherland]

Hi Nicole,

 

It might be possible to use the following to ‘un’-rasterize your raster stack into a data frame for the prediction:

 

pred.data<-stack(step1, step2, step3, sand.crop,step4, step5)  ##use rasters 'stepX' resampled to sand.crop to force extent
names(pred.data) <- c("Water.Content", "Elevation", "Habitat", "Sand.Content","NDVI", "Temperature")

## HERE

pred.dataDF <- as.data.frame(values(pred.data))

E <- predict(mAbunDA, type="state", newdata=pred.dataDF)

 

then you could make a prediction layer something for plotting like this:

 

predRast <- rasterFromXYZ(cbind(coordinates(pred.data),E))

plot(predRast)

 

or it might not…

[Richard Chandler]

That's a better idea. Thanks Chris.