# 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))