library("gamlss.dist")
data("lice")
mPO <- gamlss(head ~ 1, data = lice, family = PO, weights = freq,
  trace = FALSE)
mNBI <- gamlss(head ~ 1, data = lice, family = NBI, weights = freq,
  trace = FALSE)
mPIG <- gamlss(head ~ 1, data = lice, family = PIG, weights = freq,
  trace = FALSE)
mSI <- gamlss(head ~ 1, data = lice, family = SICHEL, weights = freq,
  method = mixed(10, 50), trace = FALSE)
AIC(mPO, mNBI, mPIG, mSI)

summary(mSI)

mSI <- histDist(lice$head, "SICHEL", freq = lice$freq, xmax = 10,
  main = "Sichel distribution", start.from = mSI, xlim = c(0, 8.75), 
  trace = FALSE)

library("gamlss.dist")
data("CD4")
plot(cd4 ~ age, data = CD4)

op <- par(mfrow = c(3, 4), mar = par("mar") + c(0, 1, 0, 0),
  pch = "+", cex = 0.45, cex.lab = 1.6, cex.axis = 1.6)
page <- c("age^-0.5", "log(age)", "age^.5", "age")
pcd4 <- c("cd4^-0.5", "log(cd4+0.1)", "cd4^.5")
for (i in 1:3) {
  yy <- with(CD4, eval(parse(text = pcd4[i])))
  for (j in 1:4) {
      xx <- with(CD4, eval(parse(text = page[j])))
      plot(yy ~ xx, xlab = page[j], ylab = pcd4[i])
  }
}
par(op)

con <- gamlss.control(trace = FALSE)
m1 <- gamlss(cd4 ~ age, sigma.fo = ~1, data = CD4, control = con)
m2 <- gamlss(cd4 ~ poly(age, 2), sigma.fo = ~1, data = CD4, control = con)
m3 <- gamlss(cd4 ~ poly(age, 3), sigma.fo = ~1, data = CD4, control = con)
m4 <- gamlss(cd4 ~ poly(age, 4), sigma.fo = ~1, data = CD4, control = con)
m5 <- gamlss(cd4 ~ poly(age, 5), sigma.fo = ~1, data = CD4, control = con)
m6 <- gamlss(cd4 ~ poly(age, 6), sigma.fo = ~1, data = CD4, control = con)
m7 <- gamlss(cd4 ~ poly(age, 7), sigma.fo = ~1, data = CD4, control = con)
m8 <- gamlss(cd4 ~ poly(age, 8), sigma.fo = ~1, data = CD4, control = con)

GAIC(m1, m2, m3, m4, m5, m7, m8)

GAIC(m1, m2, m3, m4, m5, m7, m8, k = log(length(CD4$age)))

plot(cd4 ~ age, data = CD4)
lines(CD4$age[order(CD4$age)], fitted(m7)[order(CD4$age)], col = "red")

m1f <- gamlss(cd4 ~ fp(age, 1), sigma.fo = ~1, data = CD4, control = con)
m2f <- gamlss(cd4 ~ fp(age, 2), sigma.fo = ~1, data = CD4, control = con)
m3f <- gamlss(cd4 ~ fp(age, 3), sigma.fo = ~1, data = CD4, control = con)
GAIC(m1f, m2f, m3f)

GAIC(m1f, m2f, m3f, k = log(length(CD4$age)))

m3f

m3f$mu.coefSmo

plot(cd4 ~ age, data = CD4)
lines(CD4$age[order(CD4$age)], fitted(m1f)[order(CD4$age)], col = "blue")
lines(CD4$age[order(CD4$age)], fitted(m2f)[order(CD4$age)], col = "green")
lines(CD4$age[order(CD4$age)], fitted(m3f)[order(CD4$age)], col = "red")

m2b <- gamlss(cd4 ~ bs(age), data = CD4, trace = FALSE)
m3b <- gamlss(cd4 ~ bs(age, df = 3), data = CD4, trace = FALSE)
m4b <- gamlss(cd4 ~ bs(age, df = 4), data = CD4, trace = FALSE)
m5b <- gamlss(cd4 ~ bs(age, df = 5), data = CD4, trace = FALSE)
m6b <- gamlss(cd4 ~ bs(age, df = 6), data = CD4, trace = FALSE)
m7b <- gamlss(cd4 ~ bs(age, df = 7), data = CD4, trace = FALSE)
m8b <- gamlss(cd4 ~ bs(age, df = 8), data = CD4, trace = FALSE)
GAIC(m2b, m3b, m4b, m5b, m6b, m7b, m8b)

GAIC(m2b, m3b, m4b, m5b, m6b, m7b, m8b, k = log(length(CD4$age)))

fn <- function(p) AIC(gamlss(cd4 ~ cs(age, df = p[1]), data = CD4,
  trace = FALSE), k = 2)
opAIC <- optim(par = c(3), fn, method = "L-BFGS-B", lower = c(1),
  upper = c(15))
fn <- function(p) AIC(gamlss(cd4 ~ cs(age, df = p[1]), data = CD4,
  trace = FALSE), k = log(length(CD4$age)))
opSBC <- optim(par = c(3), fn, method = "L-BFGS-B", lower = c(1),
  upper = c(15))
opAIC$par

opSBC$par

maic <- gamlss(cd4 ~ cs(age, df = 10.85), data = CD4, trace = FALSE)
msbc <- gamlss(cd4 ~ cs(age, df = 1.85), data = CD4, trace = FALSE)

m1 <- gamlss(cd4 ~ cs(age, df = 10), sigma.fo = ~cs(age, df = 2),
  data = CD4, trace = FALSE)
fn <- function(p) AIC(gamlss(cd4 ~ cs(age, df = p[1]), sigma.fo = ~cs(age,
  p[2]), data = CD4, trace = FALSE, start.from = m1), k = 2)
opAIC <- optim(par = c(8, 3), fn, method = "L-BFGS-B", lower = c(1,
  1), upper = c(15, 15))
opAIC$par

m42 <- gamlss(cd4 ~ cs(age, df = 3.72), sigma.fo = ~cs(age, df = 1.81),
  data = CD4, trace = FALSE)
m31 <- gamlss(cd4 ~ cs(age, df = 2.55), sigma.fo = ~cs(age, df = 0.93),
  data = CD4, trace = FALSE)

set.seed(1234)
rSample6040 <- sample(2, length(CD4$cd4), replace = T, prob = c(0.6,
  0.4))
fn <- function(p) VGD(cd4 ~ cs(age, df = p[1]), sigma.fo = ~cs(age,
  df = p[2]), data = CD4, rand = rSample6040)
op <- optim(par = c(3, 1), fn, method = "L-BFGS-B", lower = c(1,
  1), upper = c(10, 10))
op$par

library("gamlss.dist")
m42TF <- update(m42, family = TF)
m42PE <- update(m42, family = PE)
m42SEP3 <- update(m42, family = SEP3, method = mixed(30, 100))
m42SHASH <- update(m42, family = SHASH, method = mixed(20, 100))
GAIC(m42, m42TF, m42PE, m42SEP3, m42SHASH)

GAIC(m42, m42TF, m42PE, m42SEP3, m42SHASH, k = log(length(CD4$age)))

library"(gamlss.dist")
data("LGAclaims")
with(LGAclaims, plot(data.frame(Claims, L_Popdensity, L_KI, L_Accidents,
  L_Population)))

m0 <- gamlss(Claims ~ factor(SD) + L_Popdensity + L_KI + L_Accidents +
  L_Population, data = LGAclaims, family = PO)

m1 <- gamlss(Claims ~ factor(SD) + L_Popdensity + L_KI + L_Accidents +
  L_Population, data = LGAclaims, family = NBI)

deviance(m0) - deviance(m1)

library("MASS")
mD <- dropterm(m1, test = "Chisq")
mD

mA <- addterm(m1, scope = ~(factor(SD) + L_Popdensity + L_KI +
  L_Accidents + L_Population)^2, test = "Chisq")
mA

gs <- gamlss.scope(model.frame(Claims ~ factor(SD) + L_Popdensity +
  L_KI + L_Accidents + L_Population, data = LGAclaims))
gs

m2 <- stepGAIC.CH(m1, scope = gs, k = 2)

m2$anova

formula(m2, "mu")

op <- par(mfrow = c(3, 2))
term.plot(m2, se = T, partial = T)
par(op)

m11 <- stepGAIC.VR(m2, scope = ~L_Popdensity + L_KI + L_Accidents +
  L_Population, what = "sigma", k = 2)

library("gamlss")
data("db")
mod1 <- quote(gamlss(head~cs(nage, df = p[1]), sigma.fo = ~cs(nage, p[2]),
  nu.fo = ~cs(nage, p[3], c.spar=c(-1.5,2.5)),
  tau.fo = ~cs(nage, p[4], c.spar=c(-1.5,2.5)),
  data = db, family = BCT,
  control = gamlss.control(trace=FALSE)))
op <- find.hyper(model=mod1, other=quote(nage<-age^p[5]),
  par = c(10,2,2,2,0.25),
  lower = c(0.1,0.1,0.1,0.1,0.001),
  steps = c(0.1,0.1,0.1,0.1,0.2),
  factr = 2e9, parscale=c(1,1,1,1,0.035), penalty=2 )

library("gamlss")
data("db")
set.seed(101)
rand <- sample(2, length(db$head),replace=T, prob=c(0.6, 0.4))
table(rand)/length(db$head)

dbp <- db
dbp$agepower <- db$age^0.33
mBCT<- gamlss(head~cs(agepower,df=10.3), sigma.fo=~cs(agepower,df=3.7),
  nu.fo=~agepower, tau.fo=~agepower, family=BCT, data=dbp)

mu.s <- fitted(mBCT, "mu")[rand == 1]
sigma.s <- fitted(mBCT, "sigma")[rand == 1]
nu.s <- fitted(mBCT, "nu")[rand == 1]
tau.s <- fitted(mBCT, "tau")[rand ==1 ]
fnBCT <- function(p) {
  db$agepower <- db$age^p[1]
  vgd <- VGD(head~cs(agepower,df=p[2],c.spar=list(-1.5, 2.5)),
      sigma.fo=~cs(agepower,df=p[3],c.spar=list(-1.5, 2.5)),
      nu.fo=~agepower,  tau.fo=~agepower,
      family=BCT, data=db, rand=rand, mu.start=mu.s,
      sigma.start=sigma.s, nu.start=nu.s, tau.start=tau.s)
  cat("p=", p, " and vgd=", vgd, "\n")
  vgd }
op <- optim(par=c(.33, 12, 5.7), fnBCT, method="L-BFGS-B",
  lower=c(.01,1,1), upper=c(.5, 20, 20))

db$agepower <- db$age^0.28
mBCT <- gamlss(head~cs(agepower,df=13.77), sigma.fo=~cs(agepower,df=6.05),
  nu.fo=~agepower, tau.fo=~agepower, family=BCT, data=db)

newpar <- par(mfrow = c(2, 2), mar = par("mar") + c(0, 1, 0,0),
  col.axis = "blue4", col = "blue4", col.main = "blue4",
  col.lab = "blue4", pch = "+", cex = 0.45, cex.lab = 1.2,
  cex.axis = 1, cex.main = 1.2)
plot(mBCT, xvar = db$age, par = newpar)

par(newpar)

a <- wp(mBCT, xvar = db$age, n.inter = 20, ylim.worm = 0.6, cex = 0.3,
  pch = 20)
a

Q.stats(mBCT, xvar = db$age, n.inter = 20)

centiles.split(mBCT, xvar = db$age, c(2.5), ylab = "HEAD", xlab = "AGE",
bg= "transparent")