Error message: predictDsurface

[James Russell]

Hi,

I have been implementing a muli-session unconditional likelihood model on a fixed trapping grid (~20 m spacing) with a user defined (5 m spacing) habitat mask (3 level factor) and session covariates on Density (there is also 2 group factors, age and sex, for g(0) but I don't think they are germane to the problem).

I'd like to use predictDsurface() to get habitat-specific total estimates of density (i.e. summed across all group combinations) for each session, but on implementing predictDsurface on the fitted model I get the error:

Error in pmax(temp, 0) : cannot mix 0-length vectors with others

which I am unable to interpret

Alternatively I could try the (non-recommended) predict.secr() I suppose?

Thanks for any advice,
James

[murray.efford]

I think for your purpose the predict method is more appropriate than predictDsurface - that is for evaluating the model at every pixel across a habitat mask, whereas you are really interested only in the three levels of the habitat factor. I suspect predictDsurface (which generates the error message) is tripping up on the groups, so it's a prompt for me to actually test that part of the code, but let's bypass predictDsurface for now.

Your challenge is to construct a dataframe 'newdata' that has the requisite group, habitat and session predictors (if the default is not adequate), and then to sum the resulting predicted habitat-specific densities across groups.

I tried this with the ovenCH dataset with these results

Code: Select all

ovenmask <- make.mask(traps(ovenCH), buffer=400, spacing=15)
testfit <- secr.fit(ovenCH, mask = ovenmask, D ~ session + g, groups='Sex')
## squash the output of predict into a single dataframe
pred <- do.call(rbind, predict(testfit))
## extract the density estimates (every third row)
## two groups, so two columns
sumD <- matrix(pred$estimate[seq(1,30,3)], ncol = 2)
dimnames(sumD) <- list(2005:2009, c('F','M'))
sumD
apply(sumD,1,sum)

             F         M
2005 0.4019269 0.5180388
2006 0.4207779 0.5423356
2007 0.4972829 0.6409419
2008 0.3633985 0.4683800
2009 0.3060198 0.3944253

apply(sumD,1,sum)
     2005      2006      2007      2008      2009
0.9199657 0.9631135 1.1382248 0.8317785 0.7004452

sumDse <- matrix(pred$SE.estimate[seq(1,30,3)], ncol = 2)
dimnames(sumDse) <- list(2005:2009, c('F','M'))

apply(sumDse,1, function(x) sqrt(sum(x^2)))
    2005      2006      2007      2008      2009
0.1790869 0.1874864 0.2215747 0.1619198 0.1363535


The SE might better be combined on the link (default: log) scale, but that's a hassle. An alternative might be to change the link function for D to 'identity' in your call to secr.fit. You might want to populate a 3-D array with the predicted densities, rather than a matrix, as you have habitat, session and group as margins. Also, it is probably possible to construct the model so that it estimates the overall density directly using appropriate contrasts, but I haven't tried.

Hope this is somewhat useful, if not a complete solution.
Murray

[James Russell]

Thanks Murray. After a bit of newdata and array cobbling I got the correct density estimates for habitats. Along the way I realised I have mis-applied the groups in the secr.fit, so will need to correct the original models. This may also explain why the Dsurface was not predicting, and if it does I'll let you know.

[James Russell]

Hi,

An update with a small problem.

I created a nice new dataset for the prediction.secr function that looked like so

Code: Select all

> bagaud.newdata
    session   g seasons years hab
1         1 F.A       e     E   A
2         1 F.J       e     E   A
3         1 M.A       e     E   A
4         1 M.J       e     E   A
...

where ALL columns are factors, the first for session, the second for group, third and fourth for session covariates and the fifth for habitat (note initially I did not treat them as factors and I also annotated groups with g e.g. gF.A but neither of these appear to be the source of the problem I am having and both worked fine).

When I use this dataframe to make predictions on a secr fit model = list(D~hab+seasons+years+g) and then extract an array of density estimates using predict.secr I get correct values, which are about the same as those from the derived function (which only gives averages for habitat).

However we would also like to predict density (by habitat) from a session model such as model = list(D~hab+session+g). When we do this the estimates from predict.secr using the same newdata file differ from those using the derived function or predict directly on the model fit (without the new data). In this case the latter are correct.

I suspect the error is due to an implementation of my newdata session variable in the predict.secr function (since predict.secr with the same new data works fine for session covariates), however it is not clear to me what error I may have made. I have tried treating the session variable as both a numerical and categorical variable so I don't think this is the issue. It may be something trivial I have overlooked.

Any advice appreciated.

James

[murray.efford]

There's not much I can say about this. Parts of what you say don't make sense to me e.g. there's no way derived() (Horvitz-Thompson-like estimates) can reproduce estimates from a 'full-likelihood' model in which density is predicted from covariates. Predicting from models with factor covariates is tricky because you need to ensure the factors in 'newdata' have exactly the levels of the original covariates, in the same order.
Murray

[James Russell]

In essence my problem is, on the surface, one of imlpementing newdata incorrectly.

Specifically my problem is the implementation of my newdata session variable in predict.secr. When trying to predict with (my) newdata for a model including session the estimates of density are not correct (however predict.secr works fine if I don't specify my new dataframe, or when I use my new dataframe without session).

I can't reproduce the session error in the ovenbird or skink data, when using session or session+habclass respectively, so it must be something peculiar to my data, however I can't (for the life of me) see what it might be since the session variables in all the new dataframes are apparently identical.

I shall surrender to this one.

And yes I was implementing derived (inappropriately) on fulllikelihood models. Having worked with only conditional likelihood models in the past I didn't realise. It appeared to give ball park sensible output though, but perhaps if its a no-no there could be a stop() on it when called on full likelihood models.

Thanks.

[murray.efford]

The function 'secr.make.newdata' is used internally to create the default newdata. It also happens to be exported (i.e. revealed to users) - see ?secr.make.newdata. This is one of those rare instances when it may be worth calling it, just to compare your 'newdata' with the default. (e.g. str(secr.make.newdata(myfittedmodel)) vs str(mynewdata)).
Murray

[James Russell]

Very wise of you to include that function externally. I had wondered if I could get in to something like that to look at the guts.

Indeed there were subtleties in the newdata I couldn't isolate

Correct using secr.make.newdata

Code: Select all

> str(bagaud.newdata)
'data.frame':   144 obs. of  5 variables:
 $ session: Factor w/ 12 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ g      : Factor w/ 4 levels "F.A","F.J","M.A",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ seasons: Factor w/ 4 levels "a","e","h","p": 2 2 2 1 3 3 4 2 2 2 ...
 $ years  : Factor w/ 2 levels "E","N": 1 1 1 1 1 1 2 2 2 2 ...
 $ hab    : Factor w/ 3 levels "A","G","I": 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "out.attrs")=List of 2
  ..$ dim     : Named int  12 4
  .. ..- attr(*, "names")= chr  "session" "g"
  ..$ dimnames:List of 2
  .. ..$ session: chr  "session=1" "session=2" "session=3" "session=4" ...
  .. ..$ g      : chr  "g=F.A" "g=F.J" "g=M.A" "g=M.J"

Incorrect using as.data.frame

Code: Select all

> str(bagaud.newdata)
'data.frame':   144 obs. of  5 variables:
 $ session: Factor w/ 12 levels "1","2","3","4",..: 1 5 6 7 8 9 10 11 12 2 ...
 $ g      : Factor w/ 4 levels "F.A","F.J","M.A",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ seasons: Factor w/ 4 levels "a","e","h","p": 2 2 2 1 3 3 4 2 2 2 ...
 $ years  : Factor w/ 2 levels "E","N": 1 1 1 1 1 1 2 2 2 2 ...
 $ hab    : Factor w/ 3 levels "A","G","I": 1 1 1 1 1 1 1 1 1 1 ...

However secr.make.newdata didn't extract the habitat mask covariate levels but I just used rbind() to replicate and then alter the levels manually.

It worked and I finally have the estimates. Thanks for your patience.

[murray.efford]

It's great you were able to sort this out. Factors can be tricky, as I said. Even if I hadn't exported secr.make.newdata in the package namespace, it would be possible to see the default newdata by using the 'savenew' argument of predict.secr.

Code: Select all

>library(secr)
> temp <- predict(ovenbird.model.1, savenew = TRUE)
> attr(temp, 'newdata')
  session
1    2005
2    2006
3    2007
4    2008
5    2009


In principle you can find any function in the package source code (downloadable from http://cran.r-project.org/web/packages/secr/index.html), but that is definitely not for the faint-hearted.

Note also that the code for methods such as predict.secr can be displayed thus:

Code: Select all

> predict.secr
Error: object 'predict.secr' not found
## whereas...
> getAnywhere(predict.secr)
A single object matching ‘predict.secr’ was found
It was found in the following places
  registered S3 method for predict from namespace secr
  namespace:secr
with value

function (object, newdata = NULL, se.fit = TRUE, alpha = 0.05,
    savenew = FALSE, scaled = FALSE, ...)
{
    if (is.null(object$fit)) {
        warning("empty (NULL) object")
        return(NULL)
    }
 etc.


Murray