occu stepwise function

[Richard Schuster]

Hi all,

As there might be some interest in a stepwise procedure for the occu
function I am putting up this post to let you know that I have written
a function for this and would be happy to distribute it to anyone who
wants to use it. Please get in touch with me and I will send you the
code.

Regards,
Richard

[Tim Jessen]

Hi Richard,

My first time asking an R question to a Google group so if formatting isn't exactly as asked for, I do apologize in advance. I have successfully run all candidate models and a committee member has asked me to perform a backward elimination regression using p-values. I've seen several posts of yours on this topic, copied your code, and tried best to fit my data to it. However, I currently get an error message that states:

Error in f.AICc.occu.sig(start.model = start.model, blocks = var, max.iter = 30,  : 

  could not find function "AICc"

Would you be able to offer any help on this issue. My question is which variables (if any) influence probability of detection of this particular species at this site. Thank you so much in advance for your time and any help you may be able to provide. I have attached the dataset and my R code.

Again, thank you for any assistance you may be able to provide!

Tim

[Jeffrey Royle]

hi Tim,

 I didn't go back and read the previous posts on this topic but I think the issue is as simple as this: there is no AICc function in unmarked.

but I'm guessing this is an AICcmodavg package function so you should load that package and see what happens.

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/d/optout.

[Tim Jessen]

Thank you Andy for your help! The AICcmodavg package definitely solved the first problem. However, now, I have a second issue. It seems again like an easy fix, I just don't know how to fix it. The error message I now get is:

Error in blocks[kk] : object of type 'closure' is not subsettable

I continue to work with the code and data set I attached on my first post. If something else is needed, please ask and I'll get that to you ASAP. Thank you again for taking the time to help!

Tim

--
You received this message because you are subscribed to a topic in the Google Groups "unmarked" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/unmarked/te9mRgs1c0Q/unsubscribe.
To unsubscribe from this group and all its topics, send an email to unmarked+u...@googlegroups.com.

[Kery Marc]

Dear Tim,

your new error is a typical R error message. Bet it has nothing to do with unmarked. Go check how "blocks" look like for instance, by printing the object or doing "str(blocks)".

Best regards  --- Marc

---

From: unma...@googlegroups.com [unma...@googlegroups.com] on behalf of Tim Jessen [tim.je...@gmail.com]
Sent: 21 May 2016 03:08
To: unma...@googlegroups.com
Subject: Re: [unmarked] Re: occu stepwise function

[Dan Linden]

Tim, a quick look at your code indicates that you are assigning var to blocks in your function call:

CLCHr<- f.AICc.occu.sig(start.model=start.model, blocks=var,max.iter=30, AICcut=1)

This is what causes the error since it doesn't appear that you've defined var anywhere, which already exists as a function.  When you try to subset a function, you get the error about the closure object not being subsettable.  So, clearly you need a custom object assigned to blocks that represents something (variable names?) and it should probably be given a name that doesn't match a function.

[Tim Jessen]

Hi everyone,

Thank you all again for your help. I was able to get... something to work. Sorta. I still get error messages. I will paste the print out below. I have essentially three questions:

1. Error in model.frame.default(stateformula, siteCovs, na.action = NULL) : object is not a matrix

I received 9 of these error messages. I also ran a total of 9 models. I'm assuming it has to do something withy my site covariate data not being in a matrix. What exactly is this telling me and what is a way I can correct this?

2. Warning message: In fitList(fits = modlst) : Your list was unnamed, so model names were added as c('1','2',...)

Here I believe I understand that I did not name my models and so they were assigned 1, 2, 3... by R. Where can I go to address this issue. 

Finally, 3... did it work and what does this mean? Perhaps someone has an example of what the output is supposed to look like? The parameter I am most interested in is detection probability. All sites were occupied. I want to know how well my sampling method detects this species and what variables may/may not influence detection probability. Pasted below is my code and statistical output.

Thank you all again so much for your help!

Tim

data <- csvToUMF("2010_middens_mgrs.csv", type="unmarkedFrameOccu")

fm0 <- occu(~1 ~1, data)

fm1 <- occu(~days ~1, data)

fm2 <- occu(~mean_temp ~1, data)

fm3 <- occu(~min_temp ~1, data)

fm4 <- occu(~humidity ~1, data)

fm5 <- occu(~precip ~1, data)

fm6 <- occu(~precip*mean_temp ~1, data)

fm7 <- occu(~precip*min_temp ~1, data)

fm8 <- occu(~precip+min_temp ~1, data)

fm9 <- occu(~precip*days ~1, data)

start.model <- fm0

det.var <- as.vector(c("fm0", "fm1", "fm2", "fm3", "fm4", "fm5", "fm6", "fm7",

                       "fm8", "fm9"))

f.AICc.occu.sig <- function(start.model, blocks = det.var, max.iter=NULL,

                            detocc = 2, AICcut = 1, print.log=TRUE){

  modlst <- c(start.model)

  x <- 2

  if (detocc == 2) coeff <- length(start.model@estimates@estimates$state@estimates)

  else coeff <- length(start.model@estimates@estimates$det@estimates)

  

  best <- FALSE

  model.basis <- start.model

  keep <- list(start.model)

  

  AICmin <- AIClst <- AICc(start.model)

  

  cutoff <- 20

  

  if (is.null(max.iter) | max.iter > length(blocks)) max.iter <- length(blocks)

  

  for(ii in 1:max.iter){

    models <- list()

    coeff <- coeff + 1

    cnt <- 1

    

    for (jj in 1:length(keep)) {

      

      for(kk in 1:length(blocks)) {

        

        if (detocc == 2) form <- as.formula(paste("~. ~. + ", blocks[kk]))

        else form <- as.formula(paste("~. + ", blocks[kk], "~. "))

        

        if (class(dummy <- try(update(keep[[jj]], form))) == "unmarkedFitOccu") {

          flag <- 0

          

          if (detocc == 2) {

            

            for(dd in 1:length(sqrt(diag (vcov(dummy@estimates@estimates$state))))) {

              

              if(diag (vcov(dummy@estimates@estimates$state))[[dd]] < 0

                 || sqrt(diag (vcov(dummy@estimates@estimates$state))[[dd]]) > cutoff)

              flag <- 1

              break

            }

          }

          else {

            

            for(dd in 1:length(sqrt(diag (vcov(dummy@estimates@estimates$det))))) {

              

              if(diag (vcov(dummy@estimates@estimates$det))[[dd]] < 0

                 || sqrt(diag (vcov(dummy@estimates@estimates$det))[[dd]]) > cutoff)

              flag <- 1

              break

            }

          }

          

          if (flag == 0) {

            for (bb in 1:length(AIClst)) {

              if (round(AICc(dummy),digits=6) == round(AIClst[bb], digits=6)) {

                flag <- 1

                break

              }

            }

          }

        }

        else {flag <- 1}

        

        if (flag == 0) {

          models[[cnt]] <- dummy

          modlst[[x]] <- models[[cnt]]

          AIClst <- c(AIClst, AICc(models[[cnt]]))

          x <- x + 1

          cnt <- cnt + 1

        }

      }

    }

    

    if(length(LL <- unlist(lapply(models, function(x) {AICc(x)}))) == 0) break

    keep <- list()

    k = 1

    cont <- 0

    

    for (mm in order(LL, decreasing=FALSE)) {

      

      if (LL[mm]< AICmin + AICcut) {

        

        if (detocc == 2) {

          

          if (length(models[[mm]]@estimates@estimates$state@estimates) == coeff) {

            keep[[k]]<- models[[mm]]

            k <- k + 1

            

            if (LL[mm]<AICmin) {

              AICmin <- LL[mm]

              cont <- 1

            }

          }

        }

        else {

          if (length(models[[mm]]@estimates@estimates$det@estimates) == coeff) {

            keep[[k]]<- models[[mm]]

            k <- k + 1

            

            if (LL[mm]<AICmin) {

              AICmin <- LL[mm]

              cont <- 1

            }

          }

        }

      }

      else break

    }

    

    if (length(keep) == 0) break

  }

  

  fitlst <- fitList(fits = modlst)

  modsel <- modSel(fitlst,nullmod=NULL)

  

  return(list(model=model.basis, modlst=modlst, fitlst=fitlst, modsel=modsel))

}

MGRSr<- f.AICc.occu.sig(start.model=start.model, blocks=det.var, max.iter=30, AICcut=1)

################### TEST OUTPUT #########################

#########################################################

MGRSr

$model

Call:

  occu(formula = ~1 ~ 1, data = data)

Occupancy:

  Estimate    SE    z P(>|z|)

1.88 0.639 2.94 0.00331

Detection:

  Estimate    SE    z P(>|z|)

1.65 0.372 4.43 9.3e-06

AIC: 71.09312 

$modlst

$modlst[[1]]

Call:

  occu(formula = ~1 ~ 1, data = data)

Occupancy:

  Estimate    SE    z P(>|z|)

1.88 0.639 2.94 0.00331

Detection:

  Estimate    SE    z P(>|z|)

1.65 0.372 4.43 9.3e-06

AIC: 71.09312 

$fitlst

An object of class "unmarkedFitList"

Slot "fits":

  $`1`

Call:

  occu(formula = ~1 ~ 1, data = data)

Occupancy:

  Estimate    SE    z P(>|z|)

1.88 0.639 2.94 0.00331

Detection:

  Estimate    SE    z P(>|z|)

1.65 0.372 4.43 9.3e-06

AIC: 71.09312 

$modsel

nPars   AIC delta AICwt cumltvWt

1     2 71.09  0.00  1.00     1.00

#####################################################

###################### END ##########################

[Dylan Rhea-Fournier]

Hi Richard-
 I am interested in a stepwise procedure for the occu function if that is still available.
Thanks in advance!

cheers
Dylan

On Monday, November 29, 2010 at 9:20:08 AM UTC-7, Richard Schuster wrote:

[Richard Schuster]

Hi Dylan,

Attached the current stepwise procedure code (f.AIC_cut.occu.sig.used.15.10.23.R) and a worked example in BRCR_eBird.R.

Have a look and see if this works for you. If you need any help getting stared or adopting the code to your needs (apologies for my rather terrible annotations or lack thereof) just let me know.

Cheers,
Richard

--

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/d/optout.

-- 
Richard Schuster

Liber Ero Postdoctoral Fellow
Geomatics and Landscape Ecology Laboratory
Department of Biology
Carleton University
Email: ma...@richard-schuster.com
http://richard-schuster.com/
Twitter: @RicSchuster
ph: 250-635-2321