# read in data

age<-scan()
  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73
  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92
  93  94  95  96  97  98  99 100 101 102 103 104

popn<-scan()
    84   418  1066  2483  3721  5460  6231  8061  9487 10770 24267 26791
 29174 28476 25840 23916 21412 20116 18876 17461 15012 11871 10002  8949
  7751  6140  4718  3791  2806  2240  1715  1388   898   578   510   430
   362   291   232   196   147   100   161    11    10     8     5     4
     2     2

deaths<-scan()
   1   2  10  21  35  62  50  55  88 132 267 300 432 491 422 475 413 480 537
 566 581 464 461 433 515 374 348 304 249 167 192 171 126  86  97  93  75  84
  31  75  29  25  20   5   3   2   0   2   0   1

# define function nplogit

nplogit<-function(y,n,t,df,tol=1e-006)
{
        gt<-log((y+0.5)/(n-y+0.5))
        repeat{
                p<-1/(1+exp(-gt))
                gtnew<-smooth.spline(t,gt+(y-n*p)/(n*p*(1-p)),
                        n*p*(1-p),df=df)$y
                if(sum((gt-gtnew)^2)<tol) break
                gt<-gtnew
                }
        p
}

# plot data, call nplogit and plot fitted rurves

plot(age,deaths/popn)

library(modreg)
for(df in c(5,10,15,20))
{p<-nplogit(deaths,popn,age,df)
lines(age,p)
}

# similar analysis using built-in function gam

plot(age,deaths/popn)

library(mgcv)
deathpopn<-cbind(deaths,popn-deaths)
fit<-gam(deathpopn~s(age),binomial)
lines(age,fitted(fit))