# load data
data(faithful)

# group into 30 bins
jk <- hist(faithful$erup,breaks=seq(1.5,5.25,length=31),
          freq=TRUE,plot=FALSE)
x <- jk$mids
y <- jk$counts
# fit the model and evaluate
faithful.fit <- gssanova(y~x,family="poisson",offset=rep(log(60/8),30))
xx <- seq(1.5,5.25,length=101)
est <- predict(faithful.fit,data.frame(x=xx,offset=0),se=TRUE)
# plot the fit
plot(x,y*8/60,ylim=c(0,5.44),
     xlab="eruption duration",ylab="intensity (per sec)")
lines(xx,exp(est$fit),col=4)
lines(xx,exp(est$fit+1.96*est$se),col=5)
lines(xx,exp(est$fit-1.96*est$se),col=5)

# group into 60 bins
jk <- hist(faithful$erup,breaks=seq(1.5,5.25,length=61),
          freq=TRUE,plot=FALSE)
x <- jk$mids
y <- jk$counts
# fit the model and evaluate
faithful.fit <- gssanova(y~x,family="poisson",offset=rep(log(60/16),60))
xx <- seq(1.5,5.25,length=101)
est <- predict(faithful.fit,data.frame(x=xx,offset=0),se=TRUE)
# plot the fit
plot(x,y*16/60,ylim=c(0,5.44),
     xlab="eruption duration",ylab="intensity (per sec)")
lines(xx,exp(est$fit),col=4)
lines(xx,exp(est$fit+1.96*est$se),col=5)
lines(xx,exp(est$fit-1.96*est$se),col=5)