# Example: binomial data with lots of trials per observation
set.seed(1234)
sizevec = rep(100, length=(nn <- 200))
mydat = data.frame(x = sort(runif(nn)))
mydat = transform(mydat, prob=logit(-0+2.5*x+x^2, inverse=TRUE))
mydat = transform(mydat, y = rbinom(nn, size=sizevec, prob=prob))
mydat = transform(mydat, y = y / sizevec) # Convert to proportions
(fit = vgam(y ~ s(x, df=3), fam=amlbinomial(w=c(0.01,0.2,1,5,60)),
data=mydat, weight=sizevec, trace=TRUE))
fit@extra
par(mfrow=c(1,2))
# Quantile plot
with(mydat, plot(x, jitter(y), col="blue", las=1, main=
paste(paste(round(fit@extra$percentile, dig=1), collapse=", "),
"percentile-expectile curves")))
with(mydat, matlines(x, fitted(fit), lwd=2, col="blue", lty=1))
# Compare the fitted expectiles with the quantiles
with(mydat, plot(x, jitter(y), col="blue", las=1, main=
paste(paste(round(fit@extra$percentile, dig=1), collapse=", "),
"percentile curves are red")))
with(mydat, matlines(x, fitted(fit), lwd=2, col="blue", lty=1))
for(ii in fit@extra$percentile)
with(mydat, matlines(x, qbinom(p=ii/100, size=sizevec, prob=prob)/sizevec,
col="red", lwd=2, lty=1))
Run the code above in your browser using DataLab