www.phidot.org

I have two questions concerning AICc in RMark.

1) What could cause the AICc in the individual model results to differ from the AICc in the model.table?

I am running a multi-scale occupancy analysis. When I compare the AICc values reported in the model table

Code: Select all

 results$model.table$AICc

to the values in the individual model results

Code: Select all

results$Psi.1.Theta.1.p.1$results$AICc

I get some major discrepancies in approximately 75% of my models. Sometimes with the AICc differing by as much as 10 AICc. I tried running it with adjust=FALSE and adjust=TRUE.

2) What could cause extreme (20 AICc or more) differences between identical model runs?
If I immediately rerun models, I often get wildly different reported AICc in the results$model$results$AICc output between runs.

I have not been able to demonstrate these drastic fluctuations in simple example code so I am posting unreproducible code below with the output commented out.

Code: Select all

#1a. Process Data
proc.data <- process.data(run.data, model = "MultScalOcc", begin.time = 1, mixtures = 12)

#1b. Create the data design
ddl.run = make.design.data(proc.data)


#1c. Build a formula list within a MARK function
do.Species=function()
{
  p.1=list(formula=~TempC + RHPer + WindMS + RainMm + RecSpce + RecLnHr + YDay)
 
  Theta.1 = list(formula=~DstRoadM + DstWatrM)
 
  Psi.1 = list(formula=~FrstConct + FrstClump + FrstCorHa + FrstNDCA + DevAreaHa + WtrAreaHa + SumTempC + SumPrecMm + CliffHA)
  Psi.2 = list(formula=~FrstConct + FrstClump + FrstCorHa + FrstNDCA + DevAreaHa + WtrAreaHa + SumTempC + SumPrecMm + CliffHA + Lat + Long)
  Psi.3 = list(formula=~FrstConct + FrstClump + FrstCorHa + FrstNDCA + DevAreaHa + WtrAreaHa + SumTempC + SumPrecMm + CliffHA + Lat * Long)
 
  cml=create.model.list("MultScalOcc")
  return(mark.wrapper(cml,data=proc.data,ddl=ddl.run,adjust=FALSE,realvcv=TRUE))
}


#1d. Run the global model

a <- do.Species()
b <- do.Species()
c <- do.Species()
d <- do.Species()
e <- do.Species()
f <- do.Species()
g <- do.Species()
h <- do.Species()
i <- do.Species()
j <- do.Species()
k <- do.Species()
l <- do.Species()


a.1 <- c(a$Psi.1.Theta.1.p.1$results$AICc, a$Psi.2.Theta.1.p.1$results$AICc, a$Psi.3.Theta.1.p.1$results$AICc)
b.1 <- c(b$Psi.1.Theta.1.p.1$results$AICc, b$Psi.2.Theta.1.p.1$results$AICc, b$Psi.3.Theta.1.p.1$results$AICc)
c.1 <- c(c$Psi.1.Theta.1.p.1$results$AICc, c$Psi.2.Theta.1.p.1$results$AICc, c$Psi.3.Theta.1.p.1$results$AICc)
d.1 <- c(d$Psi.1.Theta.1.p.1$results$AICc, d$Psi.2.Theta.1.p.1$results$AICc, d$Psi.3.Theta.1.p.1$results$AICc)
e.1 <- c(e$Psi.1.Theta.1.p.1$results$AICc, e$Psi.2.Theta.1.p.1$results$AICc, e$Psi.3.Theta.1.p.1$results$AICc)
f.1 <- c(f$Psi.1.Theta.1.p.1$results$AICc, f$Psi.2.Theta.1.p.1$results$AICc, f$Psi.3.Theta.1.p.1$results$AICc)
g.1 <- c(g$Psi.1.Theta.1.p.1$results$AICc, g$Psi.2.Theta.1.p.1$results$AICc, g$Psi.3.Theta.1.p.1$results$AICc)
h.1 <- c(h$Psi.1.Theta.1.p.1$results$AICc, h$Psi.2.Theta.1.p.1$results$AICc, h$Psi.3.Theta.1.p.1$results$AICc)
i.1 <- c(i$Psi.1.Theta.1.p.1$results$AICc, i$Psi.2.Theta.1.p.1$results$AICc, i$Psi.3.Theta.1.p.1$results$AICc)
j.1 <- c(j$Psi.1.Theta.1.p.1$results$AICc, j$Psi.2.Theta.1.p.1$results$AICc, j$Psi.3.Theta.1.p.1$results$AICc)
k.1 <- c(k$Psi.1.Theta.1.p.1$results$AICc, k$Psi.2.Theta.1.p.1$results$AICc, k$Psi.3.Theta.1.p.1$results$AICc)
l.1 <- c(l$Psi.1.Theta.1.p.1$results$AICc, l$Psi.2.Theta.1.p.1$results$AICc, l$Psi.3.Theta.1.p.1$results$AICc)


data.frame(a = a.1,
           b = b.1,
           c = c.1,
           d = d.1,
           e = e.1,
           f = f.1,
           g = g.1,
           h = h.1,
           i = i.1,
           j = j.1,
           k = k.1,
           l = l.1)

#      a        b        c        d        e        f        g        h        i        j        k        l
# 1 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418 656.8418
# 2 639.4293 639.4293 639.4293 642.6255 639.4293 639.4293 639.4293 642.6255 642.6255 639.4293 642.6255 642.6255
# 3 624.0836 648.7748 621.3204 645.4115 652.2266 626.9131 652.2266 652.2266 652.2266 652.2265 652.2266 652.2266



I'd be happy to share my data off list to see if the issues could be reproduced.

Thanks

I didn't see anything in your code that would make me suspicious. My guess is that in some of the model runs MARK is not converging but it should be deterministic with the same starting values. I suggest that you look at the outputs for the models that didn't have as high an AICc to see if they converged.

1) You didn't show that the AICc values differed between model.table and what is in results for the individual model. This doesn't make sense to me. Recognize that they are listed in model.table in AICc order and not the order in which they were run.

2) This is likely a convergence issue. Possibly your data doesn't support that many parameters. It looks like model 1 gave the same results each time, so you might want to run it and then use it for starting values (see initial argument) for the other 2 models which had issues.

Whenever, you have issues like this you should look at the MARK output file and see what it is saying.

--jeff