# generate data
set.seed(5732)
test <- function(x)
        {1e+06*(x^11*(1-x)^6)+1e4*(x^3*(1-x)^10)+.1}
x <- (0:100)/100
mu <- test(x)
p <- 1/(mu/2+1)
y <- rnbinom(x,2,p)

# fit the model with known alpha
nb.fit1 <- gssanova(cbind(y,2)~x,family="nbinomial")
# fit the model with unknown alpha
nb.fit2 <- gssanova(y~x,family="nbinomial",maxiter=70)

# evaluate and plot the known alpha fit
est1 <- predict(nb.fit1,data.frame(x=x),se=TRUE)
plot(x,y,ylab=expression(alpha*(1-p(x))/p(x)))
lines(x,2*(1-p)/p,col=6)
lines(x,2/exp(est1$fit),col=4)
lines(x,2/exp(est1$fit+est1$se),col=5)
lines(x,2/exp(est1$fit-est1$se),col=5)

# evaluate and plot the unknown alpha fit
est2 <- predict(nb.fit2,data.frame(x=x),se=TRUE)
alpha <- nb.fit2$alpha
plot(x,y,ylab=expression(alpha*(1-p(x))/p(x)))
lines(x,2*(1-p)/p,col=6)
lines(x,alpha/exp(est2$fit),col=4)
lines(x,alpha/exp(est2$fit+est2$se),col=5)
lines(x,alpha/exp(est2$fit-est2$se),col=5)