gistfile1.r

library(dplyr)

# install.packages("rgeos", type="source")
# install.packages("rgbif", type="source")
# library(rgbif)
# d <- occ_search(orderKey=797, datasetKey="a8d08280-1def-11de-be11-b8a03c50a862", year="2013", limit=2000)

library(RJSONIO)


gen.template <- "http://hyonteiset.luomus.fi/insects/json?op=%s&order=Lep&%s"
kuhmo.url <- sprintf(gen.template, "search", "&county=Kuhmo&startYear=2013&endYear=2013&select=all")
jsobj <- readLines(kuhmo.url)
Encoding(jsobj) <- "latin1"
d.orig <- plyr::ldply(fromJSON(jsobj), identity)  %>%
  filter(grepl("[0-9]+W|valorysä", method) & grepl("Marko Tähtinen|Hannu Saarenmaa", observer)) %>%
  mutate(name=paste(genus, species, sep=" "), count=as.integer(totalCount)) 
d <- d %>%
  select(name, count) %>%
  group_by(name) %>%
  summarise(count=sum(count))

plot(log(-sort(-d$count)))

m.ls <- sads::fitsad(d$count, "ls")
m.ls
sqrt(m.ls@vcov[1, 1]) # approximate posterior sd of alpha

m.mzsm <- sads::fitsad(d$count, "mzsm")
m.mzsm
sqrt(m.mzsm@vcov[1, 1]) # approximate posterior sd of theta

data.frame(n.speciments=1:10,
           real.count=table(d$count)[1:10], 
           pred.ls=nrow(d)*sads::dls(1:10, sum(d$count), m.ls@coef),
           pred.mzsm=nrow(d)*sads::dmzsm(1:10, sum(d$count), m.mzsm@coef))
           
# Per yday there's too little data to get proper estimates for theta or alpha.
# (Also the gaussian approximation for the posterior breaks up.)
d2 <- d.orig %>%
  select(name, count, startDay, startMonth) %>%
  mutate(yday=30*as.integer(startMonth)+as.integer(startDay)) %>%
  group_by(yday, name) %>%
  summarise(count=sum(count))

library(ggplot2)
d2 %>% do({ m <- sads::fitsad(.$count, "ls"); data.frame(alpha=m@coef, sd=m@vcov[1, 1])})  %>% 
  ggplot(., aes(x=yday, y=alpha, ymin=alpha-sd, ymax=alpha+sd)) + geom_line() + geom_ribbon(alpha=.2)