rbchan / unmarked

R package for hierarchical models in ecological research
https://rbchan.github.io/unmarked/
37 stars 25 forks source link

Problems with predict and pcountOpen #185

Closed cschwarz-stat-sfu-ca closed 4 years ago

cschwarz-stat-sfu-ca commented 4 years ago

I have a continuous covariate for some of the parameters, but can't get the predict() function with the newData argument to get predictions. I've installed the development version of unmarked. MWE follows.

When I run the code interactively, the predicts using the newData argument fail, eg.

pred.abund1 <- predict(fitList(fits=list(m1,m2)), type = 'lambda')#, newdata = newData, appendData=TRUE) Warning message: In fitList(fits = list(m1, m2)) : Your list was unnamed, so model names were added as c('1','2',...) head(pred.abund1) Predicted SE lower upper 1 3.856037 0.5864538 2.870791 5.179436 2 4.092597 0.8570777 2.878731 5.863871 3 3.912215 0.6325138 2.883539 5.308952 4 3.818146 0.5670570 2.854541 5.107452 5 3.471858 0.8136240 2.450233 4.995011 6 3.810993 0.5647231 2.850624 5.095408

this does not work

pred.abund2 <- predict(fitList(fits=list(m1,m2)), type = 'lambda', newdata = newData, appendData=TRUE) Error: $ operator is invalid for atomic vectors In addition: Warning message: In fitList(fits = list(m1, m2)) : Your list was unnamed, so model names were added as c('1','2',...) pred.abund2 Error: object 'pred.abund2' not found

But if I run the code in Batch mode (see attachment) the code works fine and gives a prediction for the value of the site covariate=0.

Ditto for the other predictions.

Very puzzling...

Thanks Carl Schwarz test.pCount.r.txt


illustration of problems with predict() and pCountOpen models

Problems with predict

devtools::install_github("rbchan/unmarked", dependencies = TRUE, build_vignettes = TRUE) # get latest release

library(unmarked)

sessionInfo(package="unmarked")

set.seed(3) M <- 50 T <- 5 lambda <- 4 gamma <- 1.5 omega <- 0.8 p <- 0.7 y <- N <- matrix(NA, M, T) S <- G <- matrix(NA, M, T-1) N[,1] <- rpois(M, lambda) for(t in 1:(T-1)) { S[,t] <- rbinom(M, N[,t], omega) G[,t] <- rpois(M, gamma) N[,t+1] <- S[,t] + G[,t] } y[] <- rbinom(M*T, N, p)

Prepare data

umf <- unmarkedFramePCO(y = y, numPrimary=T) summary(umf)

siteCovs(umf) <- data.frame(sitecov=rnorm(M))

Fit model and backtransform

m1 <- pcountOpen(~1, ~1, ~1, ~1, umf, K=20) # Typically, K should be higher m2 <- pcountOpen(~sitecov, ~sitecov, ~sitecov, ~sitecov, umf, K=20) # Typically, K should be higher

lam <- coef(backTransform(m1, "lambda")) # or' lam

lam <- exp(coef(m1, type="lambda")) lam

gam <- exp(coef(m1, type="gamma")) gam

om<- plogis(coef(m1, type="omega")) om

p <- plogis(coef(m1, type="det")) p

ms <- modSel(fitList(fits=list(m1,m2))) ms

ms@Full

here is where things go wrong

newData<- expand.grid(sitecov=0, stringsAsFactors=TRUE) newData str(newData)

this works

pred.abund1 <- predict(fitList(fits=list(m1,m2)), type = 'lambda')#, newdata = newData, appendData=TRUE) head(pred.abund1)

this does not work

pred.abund2 <- predict(fitList(fits=list(m1,m2)), type = 'lambda', newdata = newData, appendData=TRUE) pred.abund2

this works

pred.survival1 <- predict(fitList(fits=list(m1,m2)), type = 'omega')#, newdata = newData, appendData=TRUE) head(pred.survival1)

this does not work

pred.survival2 <- predict(fitList(fits=list(m1,m2)), type = 'omega', newdata = newData, appendData=TRUE) pred.survival2

this works

pred.detect1 <- predict(fitList(fits=list(m1,m2)), type = 'det')#, newdata = newData, appendData=TRUE) head(pred.detect1)

this does not work

pred.detect2 <- predict(fitList(fits=list(m1,m2)), type = 'det', newdata = newData, appendData=TRUE) pred.detect2

this works

pred.growth1 <- predict(fitList(fits=list(m1,m2)), type = 'gamma')#, newdata = newData, appendData=TRUE) head(pred.growth1)

this does not work

pred.growth2 <- predict(fitList(fits=list(m1,m2)), type = 'gamma', newdata = newData, appendData=TRUE) pred.growth2

kenkellner commented 4 years ago

I haven't been able to replicate the issue. Your code does error for me on the release version 1.0.0 but it runs fine (including interactively) on version 1.0.0.9008. The error you're seeing is the same one that should have been fixed. I see in the batch script output that unmarked 1.0.0.9008 is definitely being loaded, can you double check that in your interactive session you have the new version loaded?

cschwarz-stat-sfu-ca commented 4 years ago

Hmm... weird... interactive session also showed 1.0.0.9008 is loaded... but still failed. But I had to a restart of my computer and now it works just fine... go figure....

Thanks for all your help in any case...

Carl