## ----setup, include = FALSE-------------------------------
knitr::opts_chunk$set(
  tidy.opts = list(width.cutoff = 60),
  collapse = TRUE,
  comment = "#>"
)
options(width = 60)

library(cNORM)

## ----pmf-figure, fig.height = 3.6, fig.width = 7, echo = FALSE----
x <- 0:45
op <- par(mar = c(4, 4, 1, 1))
plot(x, dcmp(x, mu = 20, nu = 3), type = "l", lwd = 2, col = "#1b9e77",
     ylim = c(0, 0.17), xlab = "Raw score", ylab = "P(Y = y)", bty = "l")
lines(x, dcmp(x, mu = 20, nu = 1),   lwd = 2, col = "#7570b3")
lines(x, dcmp(x, mu = 20, nu = 0.5), lwd = 2, col = "#d95f02")
legend("topright", bty = "n", lwd = 2,
       col = c("#1b9e77", "#7570b3", "#d95f02"),
       legend = c(expression(nu == 3 ~ "(under-dispersed)"),
                  expression(nu == 1 ~ "(Poisson)"),
                  expression(nu == 0.5 ~ "(over-dispersed)")))
par(op)

## ----data-------------------------------------------------
head(speed)
range(speed$fluency)

## ----dispersion-------------------------------------------
by_year <- split(speed$fluency, floor(speed$age))
round(t(sapply(by_year, function(x)
  c(n = length(x), mean = mean(x), sd = sd(x), dispersion = var(x) / mean(x)))), 2)

## ----fig1, fig.height = 4, fig.width = 7------------------
model.open <- cnorm.cmp(age = speed$age, score = speed$fluency,
  mu_degree = 3, nu_degree = 2, # default values
)

## ----ceiling-mass-----------------------------------------
age_old <- max(speed$age)
tab <- normTable.cmp(model.open, ages = age_old, start = 0, end = 150)[[1]]
p_above <- 1 - tab$Pcum[tab$x == 75]

## ----fig-ceiling, fig.height = 4, fig.width = 7-----------
model.cmp <- cnorm.cmp(
  age = speed$age,
  score = speed$fluency,
  mu_degree = 3,
  nu_degree = 2,
  max_score = 75,
  plot = TRUE
)

## ---------------------------------------------------------
summary(model.cmp, age = speed$age, score = speed$fluency)

## ----eval = FALSE-----------------------------------------
# # Grid: degrees 1 to 4 for mu, 0 to 2 for nu (0 = constant dispersion).
# # Use min_nu = -1 to include the Poisson model (nu fixed at 1) as a candidate.
# model.best <- autoselect.cmp(
#   age = speed$age,
#   score = speed$fluency,
#   max_mu = 4,
#   max_nu = 2,
#   min_nu = 0,
#   max_score = 75
# )
# 
# # Evaluated models, sorted by BIC
# head(model.best$selection$evaluated)

## ---------------------------------------------------------
cases <- data.frame(
  age = c(7.5, 8.5, 9.5, 10.5),
  raw = c(18, 24, 30, 35)
)

cases$T_Score <- round(predict(model.cmp, age = cases$age, score = cases$raw), 1)
cases$PR <- round(predict(model.cmp, age = cases$age, score = cases$raw,
                          type = "percentile"), 1)
cases

## ---------------------------------------------------------
tables <- normTable(
  c(8.0, 9.0),
  model.cmp,
  reliability = 0.88,
  CI = 0.90
)

head(tables[["8"]], 10)

## ---------------------------------------------------------
moments <- predictMoments(model.cmp, age = seq(7, 11, by = 1), censor = FALSE)
moments$dispersion <- moments$variance / moments$mean
round(moments, 2)

## ----fig2, fig.height = 4, fig.width = 7------------------
model.taylor <- cnorm(raw = speed$fluency, age = speed$age, plot = FALSE)

compare(
  model.taylor,
  model.cmp,
  age = speed$age,
  score = speed$fluency,
  title = "Taylor Polynomial vs. Conway-Maxwell-Poisson"
)