goodness of fit test: parboot- or mb.gof.test function?

[Asuncion Semper Pascual]

Hi all,

I am running a goodness of fit test to assess model fit of a single-season occupancy model. I am running the GOF test on the global model using a) the mb.gof.test function from the AICcmodavg, and b) the parboot fucntion from the unmarked package. For the parboot function I use the fitstats function from the parboot help in R to get the sums of squares, the chi-square and the Freeman-Tukey statistic.

I should get similat results, but I actually get completely different results. Here is part of the code and the output for both functions:

a) mb.gof.test

obs.boot <- mb.gof.test(Mod_global2, nsim = 1000)

print(obs.boot, digits.vals = 4, digits.chisq = 4)

MacKenzie and Bailey goodness-of-fit for single-season occupancy model

Pearson chi-square table:

                    Cohort Observed Expected Chi-square

00000000000              0        6   6.2376     0.0091

00000000010              0        1   0.2304     2.5714

00000000100              0        2   0.2304    13.5946

00000001010              0        1   0.0275    34.4165

00100010000              0        1   0.0275    34.4165

10110100001              0        1   0.0004  2325.2988

0000000001NA             1        1   0.0531    16.8734

0000100000NA             1        1   0.0531    16.8734

000000000NANA            2        3   1.8153     0.7732

00000000NANANA           3        4   2.9522     0.3718

0000000NANANANA          4        1   3.6144     1.8911

0000010NANANANA          4        1   0.2366     2.4623

0010000NANANANA          4        2   0.2366    13.1394

1000000NANANANA          4        2   0.2366    13.1394

1111111NANANANA          4        1   0.0001 16291.7951

000000NANANANANA         5        7   6.3772     0.0608

000010NANANANANA         5        1   0.4574     0.6437

000011NANANANANA         5        1   0.1156     6.7665

000100NANANANANA         5        1   0.4574     0.6437

001000NANANANANA         5        1   0.4574     0.6437

010011NANANANANA         5        1   0.0414    22.1670

00000NANANANANANA        6        6   8.5245     0.7476

00001NANANANANANA        6        1   0.6049     0.2580

00010NANANANANANA        6        1   0.6049     0.2580

00100NANANANANANA        6        3   0.6049     9.4831

01000NANANANANANA        6        2   0.6049     3.2174

0000NANANANANANANA       7        8   6.7504     0.2313

0001NANANANANANANA       7        1   0.6115     0.2468

1100NANANANANANANA       7        1   0.2225     2.7161

1111NANANANANANANA       7        1   0.0519    17.3063

000NANANANANANANANA      8       21  19.9451     0.0558

001NANANANANANANANA      8        6   3.4332     1.9190

010NANANANANANANANA      8        4   3.4332     0.0936

100NANANANANANANANA      8        1   3.4332     1.7245

110NANANANANANANANA      8        3   1.6506     1.1032

111NANANANANANANANA      8        1   0.8035     0.0480

Chi-square statistic = 18862.82 

Number of bootstrap samples = 1000

P-value = 0.008

Quantiles of bootstrapped statistics:

      0%      25%      50%      75%     100% 

   169.2    878.2   1407.0   2466.3 466944.2 

Estimate of c-hat = 6.4347 

b) parboot

fitstats <- function(Mod_global2) {

  observed <- getY(Mod_global2@data)

  expected <- fitted(Mod_global2)

  resids <- residuals(Mod_global2)

  sse <- sum(resids^2,na.rm=TRUE)

  chisq <- sum((observed - expected)^2 / expected,na.rm=TRUE)

  freeTuke <- sum((sqrt(observed) - sqrt(expected))^2,na.rm=TRUE)

  out <- c(SSE=sse, Chisq=chisq, freemanTukey=freeTuke)

  return(out)

}

pb <- parboot(Mod_global2, fitstats, nsim=10000, report=1)

cHat_pb <- pb@t0[2] / mean(p...@t.star[,2])

pb

Call: parboot(object = Mod_global2, statistic = fitstats, nsim = 10000, report = 1)

Parametric Bootstrap Statistics:

                t0 mean(t0 - t_B) StdDev(t0 - t_B) Pr(t_B > t0)

SSE           56.8           3.70             5.65       0.2559

Chisq        544.4          76.67            59.84       0.0616

freemanTukey  80.8           5.46             7.30       0.2288

t_B quantiles:

              0% 2.5% 25% 50% 75% 97.5% 100%

SSE           30   42  49  53  57    64   76

Chisq        230  373 435 464 492   589 1887

freemanTukey  44   61  70  75  80    89  103

t0 = Original statistic compuated from data

t_B = Vector of bootstrap samples

cHat_pb

   Chisq 

1.163929

As you can see in the ouput the p-value for the Chi-square test is 0.008 for the mb.gof.test function and 0.0616 for the parboot function. For c-hat I get 6.4347 for the mb.gof.test function and 1.163929 for the parboot function.

Can anybody tell me why I get such a big difference? I have read in the mb.gof.tes R documentation website the following: "Given low expected frequencies, the chi-square statistic will deviate from the theoretical distribution and it is recommended to use a parametric bootstrap approach to obtain P-values with the parboot function of the unmarked package."  

Best,

Asuncion

[Mazerolle, Marc J.]

Hi Asuncion,

mb.gof.test( ) in the AICcmodavg package implements the MacKenzie and Bailey goodness of fit test for occupancy models (see MacKenzie, D. I., and L. L. Bailey. 2004. Assessing the fit of site-occupancy models. Journal of Agricultural, Biological, and Environmental Statistics 9:300-318.), whereas the fit statistics on the raw data you are using with parboot are computed differently.

Best,

Marc
--
____________________________________
Marc J. Mazerolle
Département des sciences du bois et de la forêt
2405 rue de la Terrasse
Université Laval
Québec, Québec G1V 0A6, Canada
Tel: (418) 656-2131 ext. 7120
Email: marc.ma...@uqat.ca

________________________________________
De : unma...@googlegroups.com <unma...@googlegroups.com> de la part de Asuncion Semper Pascual <asun...@gmail.com>
Envoyé : 28 février 2018 11:02:50
À : unmarked
Objet : [unmarked] goodness of fit test: parboot- or mb.gof.test function?

--
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<mailto:unmarked+u...@googlegroups.com>.
For more options, visit https://groups.google.com/d/optout.