Problems with predict and pcountOpen
Closed this issue · 2 comments
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.095408this 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
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?
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