## Basic R code for Turk & Lorimer (2013), MNRAS, 436:3720-3726 

require(unmarked)

## PART 1: Set-up the data for analysis
## Import data and check structure

beta.data <- read.csv("/your pathway/file name.csv", header = TRUE, sep = ",")
str(beta.data)
attach(beta.data)

logLmin <- log10(Lmin)
logGamma <- log10(New.Gamma)
phat <- pnorm(logLmin, -1.1, 0.9, lower.tail = FALSE)

logphat <- log10(phat) ## Will serve as a scaled "offset"                       

beta.data <- cbind(beta.data, logLmin, phat, logphat, logGamma)

logit <- function(x) {log(x/(1-x))}
logitphat <- logit(phat)

## Figure 1
hist(Pulsars, breaks = 35, right = FALSE, ylim = c(0, 70), col = "red")

## Pull out counts
p.data <- data.frame(Pulsars)

## Create unmarkedFrame
pulsar.umf <- unmarkedFramePCount(y = p.data,
              obsCovs = data.frame(logitphat = logitphat),                    
              siteCovs = data.frame(logGamma = logGamma,
                                    Vesc = Vesc,
                                    LV = LV,
                                    phat = phat,
                                    Metal = Metal,
                                    logphat = logphat))
summary(pulsar.umf)

## PART 2: Fit models.  For brevity, we show only ZIP and NB models.

(fm1 <- pcount(~ 0 + logitphat ~ offset(phat), pulsar.umf, mixture = "ZIP", K = 200))
(fm2 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma, pulsar.umf, mixture = "ZIP", K = 200))
(fm3 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm4 <- pcount(~ 0 + logitphat ~ offset(phat) + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm5 <- pcount(~ 0 + logitphat ~ offset(phat) + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm6 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm7 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm8 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm9 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm10 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm11 <- pcount(~ 0 + logitphat ~ offset(phat) + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm12 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm13 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm14 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm15 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm16 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))

(fm17 <- pcount(~ 0 + logitphat ~ phat, pulsar.umf, mixture = "ZIP", K = 200))
(fm18 <- pcount(~ 0 + logitphat ~ phat + logGamma, pulsar.umf, mixture = "ZIP", K = 200))
(fm19 <- pcount(~ 0 + logitphat ~ phat + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm20 <- pcount(~ 0 + logitphat ~ phat + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm21 <- pcount(~ 0 + logitphat ~ phat + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm22 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm23 <- pcount(~ 0 + logitphat ~ phat + logGamma + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm24 <- pcount(~ 0 + logitphat ~ phat + logGamma + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm25 <- pcount(~ 0 + logitphat ~ phat + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm26 <- pcount(~ 0 + logitphat ~ phat + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm27 <- pcount(~ 0 + logitphat ~ phat + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm28 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm29 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm30 <- pcount(~ 0 + logitphat ~ phat + logGamma + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm31 <- pcount(~ 0 + logitphat ~ phat + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm32 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))

(fm33 <- pcount(~ 0 + logitphat ~ offset(logphat), pulsar.umf, mixture = "ZIP", K = 200))
(fm34 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma, pulsar.umf, mixture = "ZIP", K = 200))
(fm35 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm36 <- pcount(~ 0 + logitphat ~ offset(logphat) + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm37 <- pcount(~ 0 + logitphat ~ offset(logphat) + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm38 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm39 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm40 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm41 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm42 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm43 <- pcount(~ 0 + logitphat ~ offset(logphat) + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm44 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm45 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm46 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm47 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm48 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))

(fm49 <- pcount(~ 0 + logitphat ~ logphat, pulsar.umf, mixture = "ZIP", K = 200))
(fm50 <- pcount(~ 0 + logitphat ~ logphat + logGamma, pulsar.umf, mixture = "ZIP", K = 200))
(fm51 <- pcount(~ 0 + logitphat ~ logphat + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm52 <- pcount(~ 0 + logitphat ~ logphat + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm53 <- pcount(~ 0 + logitphat ~ logphat + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm54 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm55 <- pcount(~ 0 + logitphat ~ logphat + logGamma + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm56 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm57 <- pcount(~ 0 + logitphat ~ logphat + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm58 <- pcount(~ 0 + logitphat ~ logphat + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm59 <- pcount(~ 0 + logitphat ~ logphat + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm60 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm61 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm62 <- pcount(~ 0 + logitphat ~ logphat + logGamma + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm63 <- pcount(~ 0 + logitphat ~ logphat + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm64 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))

(fm65 <- pcount(~ 0 + logitphat ~ 1, pulsar.umf, mixture = "ZIP", K = 200))
(fm66 <- pcount(~ 0 + logitphat ~ logGamma, pulsar.umf, mixture = "ZIP", K = 200))
(fm67 <- pcount(~ 0 + logitphat ~ Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm68 <- pcount(~ 0 + logitphat ~ LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm69 <- pcount(~ 0 + logitphat ~ Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm70 <- pcount(~ 0 + logitphat ~ logGamma + Vesc, pulsar.umf, mixture = "ZIP", K = 200))
(fm71 <- pcount(~ 0 + logitphat ~ logGamma + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm72 <- pcount(~ 0 + logitphat ~ logGamma + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm73 <- pcount(~ 0 + logitphat ~ Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm74 <- pcount(~ 0 + logitphat ~ Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm75 <- pcount(~ 0 + logitphat ~ LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm76 <- pcount(~ 0 + logitphat ~ logGamma + Vesc + LV, pulsar.umf, mixture = "ZIP", K = 200))
(fm77 <- pcount(~ 0 + logitphat ~ logGamma + Vesc + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm78 <- pcount(~ 0 + logitphat ~ logGamma + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm79 <- pcount(~ 0 + logitphat ~ Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))
(fm80 <- pcount(~ 0 + logitphat ~ logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "ZIP", K = 200))

## Now do NB models

(fm81 <- pcount(~ 0 + logitphat ~ offset(phat), pulsar.umf, mixture = "NB", K = 200))
(fm82 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma, pulsar.umf, mixture = "NB", K = 200))
(fm83 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm84 <- pcount(~ 0 + logitphat ~ offset(phat) + LV, pulsar.umf, mixture = "NB", K = 200))
(fm85 <- pcount(~ 0 + logitphat ~ offset(phat) + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm86 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm87 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + LV, pulsar.umf, mixture = "NB", K = 200))
(fm88 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm89 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm90 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm91 <- pcount(~ 0 + logitphat ~ offset(phat) + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm92 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm93 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm94 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm95 <- pcount(~ 0 + logitphat ~ offset(phat) + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm96 <- pcount(~ 0 + logitphat ~ offset(phat) + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))

(fm97 <- pcount(~ 0 + logitphat ~ phat, pulsar.umf, mixture = "NB", K = 200))
(fm98 <- pcount(~ 0 + logitphat ~ phat + logGamma, pulsar.umf, mixture = "NB", K = 200))
(fm99 <- pcount(~ 0 + logitphat ~ phat + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm100 <- pcount(~ 0 + logitphat ~ phat + LV, pulsar.umf, mixture = "NB", K = 200))
(fm101 <- pcount(~ 0 + logitphat ~ phat + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm102 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm103 <- pcount(~ 0 + logitphat ~ phat + logGamma + LV, pulsar.umf, mixture = "NB", K = 200))
(fm104 <- pcount(~ 0 + logitphat ~ phat + logGamma + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm105 <- pcount(~ 0 + logitphat ~ phat + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm106 <- pcount(~ 0 + logitphat ~ phat + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm107 <- pcount(~ 0 + logitphat ~ phat + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm108 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm109 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm110 <- pcount(~ 0 + logitphat ~ phat + logGamma + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm111 <- pcount(~ 0 + logitphat ~ phat + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm112 <- pcount(~ 0 + logitphat ~ phat + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))

(fm113 <- pcount(~ 0 + logitphat ~ offset(logphat), pulsar.umf, mixture = "NB", K = 200))
(fm114 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma, pulsar.umf, mixture = "NB", K = 200))
(fm115 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm116 <- pcount(~ 0 + logitphat ~ offset(logphat) + LV, pulsar.umf, mixture = "NB", K = 200))
(fm117 <- pcount(~ 0 + logitphat ~ offset(logphat) + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm118 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm119 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + LV, pulsar.umf, mixture = "NB", K = 200))
(fm120 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm121 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm122 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm123 <- pcount(~ 0 + logitphat ~ offset(logphat) + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm124 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm125 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm126 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm127 <- pcount(~ 0 + logitphat ~ offset(logphat) + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm128 <- pcount(~ 0 + logitphat ~ offset(logphat) + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))

(fm129 <- pcount(~ 0 + logitphat ~ logphat, pulsar.umf, mixture = "NB", K = 200))
(fm130 <- pcount(~ 0 + logitphat ~ logphat + logGamma, pulsar.umf, mixture = "NB", K = 200))
(fm131 <- pcount(~ 0 + logitphat ~ logphat + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm132 <- pcount(~ 0 + logitphat ~ logphat + LV, pulsar.umf, mixture = "NB", K = 200))
(fm133 <- pcount(~ 0 + logitphat ~ logphat + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm134 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm135 <- pcount(~ 0 + logitphat ~ logphat + logGamma + LV, pulsar.umf, mixture = "NB", K = 200))
(fm136 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm137 <- pcount(~ 0 + logitphat ~ logphat + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm138 <- pcount(~ 0 + logitphat ~ logphat + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm139 <- pcount(~ 0 + logitphat ~ logphat + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm140 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm141 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm142 <- pcount(~ 0 + logitphat ~ logphat + logGamma + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm143 <- pcount(~ 0 + logitphat ~ logphat + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm144 <- pcount(~ 0 + logitphat ~ logphat + logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))

(fm145 <- pcount(~ 0 + logitphat ~ 1, pulsar.umf, mixture = "NB", K = 200))
(fm146 <- pcount(~ 0 + logitphat ~ logGamma, pulsar.umf, mixture = "NB", K = 200))
(fm147 <- pcount(~ 0 + logitphat ~ Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm148 <- pcount(~ 0 + logitphat ~ LV, pulsar.umf, mixture = "NB", K = 200))
(fm149 <- pcount(~ 0 + logitphat ~ Metal, pulsar.umf, mixture = "NB", K = 200))
(fm150 <- pcount(~ 0 + logitphat ~ logGamma + Vesc, pulsar.umf, mixture = "NB", K = 200))
(fm151 <- pcount(~ 0 + logitphat ~ logGamma + LV, pulsar.umf, mixture = "NB", K = 200))
(fm152 <- pcount(~ 0 + logitphat ~ logGamma + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm153 <- pcount(~ 0 + logitphat ~ Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm154 <- pcount(~ 0 + logitphat ~ Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm155 <- pcount(~ 0 + logitphat ~ LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm156 <- pcount(~ 0 + logitphat ~ logGamma + Vesc + LV, pulsar.umf, mixture = "NB", K = 200))
(fm157 <- pcount(~ 0 + logitphat ~ logGamma + Vesc + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm158 <- pcount(~ 0 + logitphat ~ logGamma + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm159 <- pcount(~ 0 + logitphat ~ Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))
(fm160 <- pcount(~ 0 + logitphat ~ logGamma + Vesc + LV + Metal, pulsar.umf, mixture = "NB", K = 200))

fms <- fitList(fm1, fm2, fm3, fm4, fm5, fm6, fm7, fm8, fm9, 
               fm10, fm11, fm12, fm13, fm14, fm15, fm16, fm17, fm18, fm19,
               fm20, fm21, fm22, fm23, fm24, fm25, fm26, fm27, fm28, fm29,
               fm30, fm31, fm32, fm33, fm34, fm35, fm36, fm37, fm38, fm39,
               fm40, fm41, fm42, fm43, fm44, fm45, fm46, fm47, fm48, fm49,
               fm50, fm51, fm52, fm53, fm54, fm55, fm56, fm57, fm58, fm59,
               fm60, fm61, fm62, fm63, fm64, fm65, fm66, fm67, fm68, fm69,
               fm70, fm71, fm72, fm73, fm74, fm75, fm76, fm77, fm78, fm79,
               fm80, fm81, fm82, fm83, fm84, fm85, fm86, fm87, fm88, fm89,
               fm90, fm91, fm92, fm93, fm94, fm95, fm96, fm97, fm98, fm99,
               fm100, fm101, fm102, fm103, fm104, fm105, fm106, fm107, fm108, fm109,
               fm110, fm111, fm112, fm113, fm114, fm115, fm116, fm117, fm118, fm119,
               fm120, fm121, fm122, fm123, fm124, fm125, fm126, fm127, fm128, fm129,
               fm130, fm131, fm132, fm133, fm134, fm135, fm136, fm137, fm138, fm139,
               fm140, fm141, fm142, fm143, fm144, fm145, fm146, fm147, fm148, fm149,
               fm150, fm151, fm152, fm153, fm154, fm155, fm156, fm157, fm158, fm159,
               fm160)

## Rank them by AIC
(ms <- modSel(fms))

## Most parsimonious model
(fm146 <- pcount(~ 0 + logitphat ~ logGamma, pulsar.umf, mixture = "NB", K = 500))   
summary(fm146)  

## logGamma effect estimate and CI on data scale
(lc <- linearComb(fm146['state'], c(0, 1)))
exp(lc@estimate)  
exp(confint(lc))  

## PART 3: Do some analysis of the results 

## Expected abundance over range of logGamma
oldpar <- par(mar=c(5,5,4,2) + 0.1)
summary(logGamma)
newData1 <- data.frame(logGamma = seq(-3, 4, length = 94))  
E.N <- predict(fm146, type = "state", newdata = newData1, appendData = TRUE)

## Plot predictions with 95% CI to make Figure 2
plot(Predicted ~ logGamma, E.N, type = "l", ylim = c(-.1,max(E.N$Predicted)),
     xlab = "Log Base 10 of Gamma",
     ylab = expression("Estimated mean number of pulsars,"~widehat(lambda))) 
lines(lower ~ logGamma, E.N, type = "l", col = gray(0.5))
lines(upper ~ logGamma, E.N, type = "l", col = gray(0.5))
points(logGamma, rep(0.0, 94))

## Expected detection probability as function of logitphat
summary(logitphat)
newData2 <- data.frame(logitphat = seq(-7.7, 0.3, length = 94))
E.p <- predict(fm146, type = "det", newdata = newData2, appendData = TRUE)
 
## Figure 3 
plot(Predicted ~ logitphat, E.p, type = "l", ylim = c(0, 1),
     xlab = "Logit of empirical estimated detection probability", xlim = c(-8, 0),
     ylab = expression("Estimated detection probability,"~widehat(p)))
lines(lower ~ logitphat, E.p, type = "l", col = gray(0.5))
lines(upper ~ logitphat, E.p, type = "l", col = gray(0.5))
points(logitphat, rep(0.0, 94))
par(oldpar)

## Goodness-of-Fit
## Function returning three fit-statistics.
fitstats <- function(fm) {
    observed <- getY(fm@data)
    expected <- fitted(fm)
    resids <- residuals(fm)
    sse <- sum(resids^2)
    chisq <- sum((observed - expected)^2 / expected)
    freeTuke <- sum((sqrt(observed) - sqrt(expected))^2)
    out <- c(SSE = sse, Chisq = chisq, freemanTukey = freeTuke)
    return(out)
    }

(pb <- parboot(fm146, fitstats, nsim = 999, report = 1))
 
(null.nb <- pcount(~ 1 ~ 1, pulsar.umf, mixture = "NB", K = 500))
summary(null.nb) 

(null.pb <- parboot(null.nb, fitstats, nsim = 999, report = 1))

## Bootstrap a summary statistic, estimated total population size (derived) 
sum(predict(fm146, type = "state")$Predicted, na.rm = TRUE)   

Nhat <- function(fm) {
    N <- sum(predict(fm, type = "state")$Predicted, na.rm = TRUE)
    }
    
(pb.N <- parboot(fm146, Nhat, nsim = 999, report = 5))

## Figure 4
plot(pb.N)   
dens <- density(pb.N@t.star)
h <- hist(pb.N@t.star, plot = FALSE)
xlim <- range(dens$x)
ylim <- range(0, h$density, dens$y)
hist(pb.N@t.star, main = "Parametric Bootstrap Distribution of 
     1,000 Estimates of N", xlab = "N", xlim = xlim, ylim = ylim, probability = TRUE)
lines(dens, lwd = 2)

## Estimates of conditional abundance posterior distribution at each GC
## with posterior mean and 95% CI for N

(re <- ranef(fm146))
sum(bup(re))   
colSums(confint(re))   

## Figure 5
p1 <- plot(re, subset = site %in% c(44), layout=c(1, 1), xlim = c(50, 200))  
dimnames(p1)$site.c <- "Ter 5"
plot(p1)

## Lambda-hat for Table 3 
cbind(beta.data[,1], predict(fm146, type = "state"))

## Figure 6  
f <- function(x){x^{0.69}}
plot(f, 1, 500, log = "xy", ylim = c(1,100), xlab = "Gamma", ylab = "Number of pulsars") 
mine <- function(x){exp(1.5*log10(x))}
curve(mine, add = TRUE)