Nmix.gof.test vs. parboot()

[Christopher F. Hansen]

Hi unmarked Google group,

I'll start off by saying that I am a novice when it comes to using R and even more so with unmarked...so make no assumptions about my knowledge!  That being said, I have two questions:

1) What is the difference between using the Nmix.gof.test() function and the parboot() function?  For example:

pbfm19 <- Nmix.gof.test(fm19,nsim=1000,plot.hist=TRUE)  ## fm19 is my model

vs.

fitstats <- function(fm19) {
  observed <- getY(BTNW_umf)  ## BTNW_umf is my data with covariates
  expected <- fitted(fm19) ## fm19 is my model
  resids <- residuals(fm19)
  sse <- sum(resids^2) # Sums of squares
  chisq <- sum((observed - expected)^2 / expected) # Chisquared
  freeTuke <- sum((sqrt(observed) - sqrt(expected))^2) # Freeman-Tukey
  out <- c(SSE=sse, Chisq=chisq, freemanTukey=freeTuke)
  return(out)
}
pbfm19 <- parboot(fm19,fitstats, nsim = 1000, report = 1)
plot(pbfm19)

I originally thought that these were the same, however, the bootstrapped distributions are very different.  I realize that with random sampling that they would be slightly different, but with 1000 simulations Nmix.gof.test() gives a normal distribution (see attached Figure 1. Nmix.gof.test) and parboot() gives a very right skewed distribution for simulated chi-square values (see attached Figure 2. parboot()).  Chi-square test statistics for the original data are the same for both methods as would be expected, however, the goodness of fit test fails for Nmix.gof.test() and passes for parboot().  Is Nmix.gof.test() a non-parametric bootstrap?  Clarification on this would be greatly appreciated. 

2) My second question is about c-hat.  My understanding is that a c-hat value greater than one potentially indicates that the original data is over-dispersed and that it can be used as a correction factor for the variance.  For my models, c-hat values range from 0.4 - 0.6.  Is this bad?  Does this indicate that the data are under-dispersed and is there something I should/can do about that?

Thanks,
Chris

[Marc J. Mazerolle]

Hi Christopher,

Do you have any missing values in your data? Essentially
Nmix.gof.test( ) just calls parboot( ) with a Pearson chi-square
computed by Nmix.chisq( ) for each iteration, and then computes the
c-hat and formats the output. The statistics are computed after removing
NA's. Nmix.gof.test( ) uses a parametric bootstrap, so you should get
the same answer as the "manual approach" you are using. If you want, you
can send me the data and minimal R code off list to reproduce the
problem and I can look into it.

Regarding the c-hat < 1, usually you would not correct for
underdispersion by changing the c-hat value. However, low values can
also identify lack of fit and you might want to try a different model
instead.

Sincerely,

Marc




--
____________________________________
Marc J. Mazerolle
Centre d'étude de la forêt
Université du Québec en Abitibi-Témiscamingue
445 boulevard de l'Université
Rouyn-Noranda, Québec J9X 5E4, Canada
Tel: (819) 762-0971 ext. 2458
Email: marc.ma...@uqat.ca


-------- Message initial --------
De: Christopher F. Hansen <chris.f...@gmail.com>
À: unma...@googlegroups.com
Sujet: [unmarked] Nmix.gof.test vs. parboot()
Date: Fri, 6 Jun 2014 11:16:22 -0700

[Christopher F. Hansen]

Thank you Marc for your response.  From what I can tell, it appears that when running the "manual approach" (i.e., creating the function and then running parboot()) you get similar results to Nmix.gof.test() if you don't put in your model name and unmarked data frame into the created function.

For example:

fitstats <- function(fm) {
  observed <- getY(fm@data)
  expected <- fitted(fm)
  resids <- residuals(fm)

  sse <- sum(resids^2) # Sums of squares
  chisq <- sum((observed - expected)^2 / expected) # Chisquared
  freeTuke <- sum((sqrt(observed) - sqrt(expected))^2) # Freeman-Tukey
  out <- c(SSE=sse, Chisq=chisq, freemanTukey=freeTuke)
  return(out)
}

as opposed to:

fitstats <- function(fm14) { # fm14 is the model with the lowest AIC
  observed <- getY(BTNW_umf) # BTNW_umf is the unmarked frame with raw data and covariates
  expected <- fitted(fm14)
  resids <- residuals(fm14)

  sse <- sum(resids^2) # Sums of squares
  chisq <- sum((observed - expected)^2 / expected) # Chisquared
  freeTuke <- sum((sqrt(observed) - sqrt(expected))^2) # Freeman-Tukey
  out <- c(SSE=sse, Chisq=chisq, freemanTukey=freeTuke)
  return(out)
}

I've seen both ways done within this forum and in online modules/examples.  At least with my data, they produce very different results.  Which ways is correct?

Thanks,
Chris

[Marc.Ma...@uqat.ca]

Hi Chris,

If you call directly the original data set each time in the function with getY(BTNW_umf), this will create problems. You should be extracting the simulated data at each iteration with getY(fm@data). The way the goodness of fit test works is that for each iteration, you (or parboot( )) simulate a data set and run the model on the simulated data. The chi-square (or other statistic) for a given iteration is computed by comparing the observed data (i.e., the data simulated at that iteration) to the fitted values (obtained from fitting the model to the simulated data at this same iteration). By repeating the procedure many many times, you get a theoretical distribution for the test statistic. You can then compare the original test statistic obtained from the original data set to the theoretical distribution. Essentially, the P-value reflects how likely it is to observe values as large or larger than the one you observed with the original data when the model fits the data.

Best,

Marc

--
____________________________________
Marc J. Mazerolle
Centre d'étude de la forêt
Université du Québec en Abitibi-Témiscamingue
445 boulevard de l'Université
Rouyn-Noranda, Québec J9X 5E4, Canada
Tel: (819) 762-0971 ext. 2458
Email: marc.ma...@uqat.ca

---

De : unma...@googlegroups.com [unma...@googlegroups.com] de la part de Christopher F. Hansen [chris.f...@gmail.com]
Envoyé : lundi 9 juin 2014 16:56
À : unma...@googlegroups.com
Objet : [unmarked] Re: Nmix.gof.test vs. parboot()

--
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.

[Giancarlo Sadoti]

Greetings unmarkers,

This seems like a relevant-enough thread to post this (I hope), and perhaps mainly a question for Marc M. (but any help is welcome!)

I've been running Nmix.gof.test on a 30+ year model from a data set with a max count of 39. I originally used K=50, which takes 10 or so minutes to fit on a newer machine.

I was interested in the getting the fit via a gof test but had issues with Nmix.gof.test returning "Error in pcountOpen(lambdaformula = ~1, gammaformula = ~1, omegaformula = ~1,  :
specified K is too small. Try a value larger than any observation".  I tried upping K to 70 (which takes about 30 minutes per model) and used 100 sims (I know more would be ideal). 

Two days in I got the same error, and sadly, no chi-square values from any of the (perhaps 50?) simulations. Raising the value of K further is a (painfully slow) possibility, but I'd *really* like to avoid the same fate when the simulation randomly generates a value above K.

Can anyone suggest a work-around?  I can see Marc M.'s code for Nmix.gof.test() (though not for parboot()), but am unsure how to modify it so that it simply skips a simulation (or returns NA) if K > max count.

Many thanks.

Giancarlo

[Marc J. Mazerolle]

Giancarlo,

The error is returned from parboot( ) called internally from
Nmix.gof.test( ). In the original model with the actual data, did you
see changes in the estimates or log-Likelihood when you increased K?
After some value of K, the estimates should stop changing. Otherwise, it
may indicate that K = 50 is not large enough. It also depends on the
detection probability - if it is low, then K = 50 might not be large
enough.

Regarding throwing away the iterations when there are problems of
convergence, Richard would be the best to address this issue.

Best,

Marc
--
____________________________________
Marc J. Mazerolle
Centre d'étude de la forêt
Université du Québec en Abitibi-Témiscamingue
445 boulevard de l'Université
Rouyn-Noranda, Québec J9X 5E4, Canada
Tel: (819) 762-0971 ext. 2458
Email: marc.ma...@uqat.ca

-------- Message initial --------
De: Giancarlo Sadoti <gcsa...@gmail.com>
À: unma...@googlegroups.com
Sujet: [unmarked] Re: Nmix.gof.test vs. parboot()
Date: Tue, 7 Oct 2014 12:24:48 -0700

[Richard Chandler]

Hi Giancarlo,

Unfortunately, these error messages are indicating that K really is too low, and the only solution is to increase K. If you discarded the simulations where K was too low, your results wouldn't be valid. 

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 Chandler

Assistant Professor

Warnell School of Forestry and Natural Resources

University of Georgia

706-542-5815

chandlerlab.wordpress.com

[Ted Jones]

Hi Giancarlo,

I'm not familiar with the Nmix.gof.test() function, but for parboot() you  can easily run your simulations as a parallel process.  Here's a link where Chandler helped me accomplish that task:

https://groups.google.com/forum/#!searchin/unmarked/parallel/unmarked/txp1UeptjgM/kOJtMWvSxTsJ

Since you've got a new computer, I assume it is at least a quad core, and thankfully R takes full advantage of logical cores, therefore, if it is hyperthreaded (like the i7) you will have twice as many virtual cores to work with.  I've had luck running parboot simulations on n-1 cores, where n is the number of virtual, or "logical", cores (eg, the i7 is a hyperthreaded quad core CPU, with a total of 8 virtual cores, so I usually parallel process on 7 cores).  

Parallel processing will not speed up the iteration time, but at least, for a given value of K, you might get 7x as many simulations accomplished in the same amount of time.  My co-worker occasionally lends me his hyperthreaded 6-core machine, and I am so glad to get 11 returns for a single loop when the lap time is measured in days!

I'm very grateful with Chandler's help on this-- what a game changer.

Also, I believe the code for parboot() is available on Chandler's GitHub page. 

Best,

Ted

[rset...@vt.edu]

Hi,

I am hoping to get some clarity on what K actually is.

I too am getting the "specified K is too small" error when running the Nmix.gof.test.  To troubleshoot, I took the parboot function given above and am having it print the observed data within the loop.  I have K=20 and get the error on this data draw :

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]

[1,]    1    0    2    1    2    2    2    4    4     5     1     4     0     0

[2,]    4    3    1    2    4    4    2    1    2     3     1     2     4     3

[3,]    0    0    0    0    4    4    1    4    4     2     7     4     1     1

[4,]    7    3    5    1    1    1    2    1    1     1     1     1     3     0

[5,]    7    3    2    1    1    5    0    1    3     1     1     1     1     0

     [,15]

[1,]     4

[2,]     2

[3,]     0

[4,]     0

[5,]     1

I do believe I have a K larger than any of this data.  The largest value in the original dataset is 10.  How should I choose K to avoid this?  I have tried to jack up K to 100 etc, but am not able to get 1000 sims in before getting the error to get a good bootstrap of the data.

Here is the original data, perhaps with this data I shouldn't expect to get to 1000 bootstrapped draws?

  y.1 y.2 y.3 y.4 y.5 y.6 y.7 y.8 y.9 y.10 y.11 y.12 y.13 y.14 y.15     

1   2   1   1   3   0   1   0   1   4    5    3    4    3    2    1 

2   7   3   3   3   5   2   3   2   1    4    5    0    2    2    1 

3   0   1   0   2   0   0   1   3   3    1    3    1    1    3    1 

4  10   6   1   0   1   0   1   1   3    2    0    1    1    1    2

5   4   2   4   1   0   0   1   0   1    3    1    2    0    1    2

[Jeffrey Royle]

the N-mixture model likelihood is formulated by summing the binomial likelihood for the counts over all possible values of N = max(count):infinity except that we approximate infinity by K, in order to do the summation. So usually K = 50 or 100 or 200 is sufficient to obtain stable MLEs but sometimes , especially for the negative binomial model, the MLEs are very unstable even for very large values of K.

The parboot function is simulating data sets but using the K from the fitted model, and it seems that in your case one of the simulated data sets has an observed count larger than K.   This might be fixed by setting K = 200 or 400 or something but, on the other hand, the likelihood may be unstable and so the problem will persist if that is the case. 

--
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+unsubscribe@googlegroups.com.

[rset...@vt.edu]

Thank you for the reply.

We ran the pcountopen with k=20 which is 2*max(count), but when doing a gof test, I keep getting the too small error.  Transiently, it finishes, so I suppose we are just getting lucky with a stable set of draws.

For valid gof test numbers, do we need to match the k in the gof test to that used in pcountopen?

To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.