#
# This is R or SPlus code for the insurance data example
#
#
# Begin by reading in the data from a file called insurance
# The data are stored in a matrix called x
#
x_read.table("insurance",header=T)
#
# Now create a new variable called Code by subtracting 1975.5
# from the variable Year. Add a 3rd column to x to hold Code
#
x_cbind(x,Code=x$Year-1975.5)
#
# To make a graph in R or SPlus you should specify a plotting device
# In this example I save the graph in a postscript file
#
postscript("dataplot.ps")
#
# plot Code as the x variable (3rd column of x) against Cost (2nd column of x)
# label the axes (using xlab and ylab), add a title (using main), make the
# axis for Code run from -5 to 6 (xlim) and don't print the actual x axis
#
plot(x[,3],x[,2],main="Claims per policy: NB 1971-1980",xlab="Year",ylab="Cost ($)",xlim=c(-5,6),xaxt='n')
#
# now add a customized x axis with ticks where I want them (at) and labelled not
# by the value of code but by the corresponding value of Year
#
axis(1,at=c(-4.5,-2.5,-0.5,1.5,3.5,5.5),labels=c(1971,1973,1975,1977,1979,1981))
#
# Tell R or SPlus that the graph is finished
dev.off()
#
# Same graph but in portrait orientation
#
postscript("classdataplot.ps",horizontal=F)
plot(x[,3],x[,2],main="Claims per policy: NB 1971-1980",xlab="Year",ylab="Cost ($)",xlim=c(-5,6),xaxt='n')
axis(1,at=c(-4.5,-2.5,-0.5,1.5,3.5,5.5),labels=c(1971,1973,1975,1977,1979,1981))
dev.off()
#
# Now add to the plots graphs for fitted values for various models
#
postscript("dataplusfits.ps")
plot(x[,3],x[,2],main="Claims per policy: NB 1971-1980",xlab="Year",ylab="Cost ($)",xaxt='n',xlim=c(-5,7),ylim=c(40,240))
axis(1,at=c(-4.5,-2.5,-0.5,1.5,3.5,5.5),labels=c(1971,1973,1975,1977,1979,1981))
#
# Fit a linear model in R or SPlus using lm which functions like glm -- more or less
# the statement model Cost = Code in SAS becomes
# Cost ~ Code
fit1_lm(Cost ~Code, data=x)
fit2_lm(Cost ~Code+Code^2, data=x)
fit3_lm(Cost ~Code+Code^2+Code^3, data=x)
fit4_lm(Cost ~Code+Code^2+Code^3+Code^4, data=x)
fit5_lm(Cost ~Code+Code^2+Code^3+Code^4+Code^5, data=x)
#
# Plot a straight line with coefficients from the straightline fit
#
abline(fit1)
#
# for models which are not straightlines we must work harder to plot
# the fitted response. Here is a function which plots a curve
# for one of our fits -- it graphs the function between x=l and x=u
# and assumes the fit is a polynomial.
#
curve_function(l,u,npts=250,fit){
cc_coef(fit)
xx_seq(l,u,length=npts)
p_length(cc)-1
xm_outer(xx,0:p,"^")
yf_xm%*%cc
list(x=xx,y=yf)
}
# That was the end of the function
#
# Now plot the curve: lty chooses different types of dotted and dashed lines.
#
lines(curve(-5,7,,fit2),lty=2)
lines(curve(-5,7,,fit3),lty=3)
lines(curve(-5,7,,fit4),lty=4)
lines(curve(-5,7,,fit5),lty=5)
#
# add a legend. x and y control where the top left of the legend box goes on
# the plot
#
legend(x=-4.5,y=200,lty=1:5,
legend=c("Degree = 1","Degree = 2","Degree = 3","Degree = 4","Degree = 5"))
#
# add a vertical line at year=1982.25 which is code 6.75
#
abline(v=1982.25-1975.5)
dev.off()
#
# Repeat in Portrait format
#
postscript("classdataplusfits.ps",horizontal=F)
plot(x[,3],x[,2],main="Claims per policy: NB 1971-1980",xlab="Year",ylab="Cost ($)",xaxt='n',xlim=c(-5,7),ylim=c(40,240))
axis(1,at=c(-4.5,-2.5,-0.5,1.5,3.5,5.5),labels=c(1971,1973,1975,1977,1979,1981))
fit1_lm(Cost ~Code, data=x)
fit2_lm(Cost ~Code+Code^2, data=x)
fit3_lm(Cost ~Code+Code^2+Code^3, data=x)
fit4_lm(Cost ~Code+Code^2+Code^3+Code^4, data=x)
fit5_lm(Cost ~Code+Code^2+Code^3+Code^4+Code^5, data=x)
abline(fit1)
lines(curve(-5,7,,fit2),lty=2)
lines(curve(-5,7,,fit3),lty=3)
lines(curve(-5,7,,fit4),lty=4)
lines(curve(-5,7,,fit5),lty=5)
abline(v=1982.25-1975.5)
legend(x=-4.5,y=200,lty=1:5,
legend=c("Degree = 1","Degree = 2","Degree = 3","Degree = 4","Degree = 5"))
dev.off()