#########################################################################
###                                                                   ###
### This Splus program generates Figure 1 used in our paper "Modeling ###
### Daily and Subdaily Cycles in Rat Sleep Data" by Peihua Qiu, Rick  ###
### Chappell, William Obermeyer and Ruth Benca. This article will be  ###
### published in Biometrics as a consultant forum paper. All questions###
### about the program should be forwarded to:	                      ###
###                                                                   ###
###                        Peihua Qiu                                 ###
###                        School of Statistics                       ###
###                        University of Minnesota                    ###
###                        313 Ford Hall                              ###
###                        224 Church Street S.E.                     ###
###                        Minneapolis, MN 55455                      ###
###                        Office Phone: (612) 625-5819               ###
###                        E-Mail: qiu@stat.umn.edu                   ###
###                                                                   ###
### @ COPYRIGHT 1999 BY PEIHUA QIU                                    ###
###                                                                   ### 
#########################################################################

postscript("figure1.ps", width=7, height=8, horizontal=F)

par(mfrow=c(2,2),oma=c(1,0,6,1),mar=c(4,4,2,1))


### Read baseline data from data file "baseline.dat" and test data from
### data file "test.dat". Also generate "time" variable which denotes 
### the middle point of each 5-minute interval of the day. The "time"
### variable is used as an independent variable.

bhbdat <- scan("baseline.dat")
bhtdat <- scan("test.dat")
time   <- seq(2.5,by=5.0,length=288)/60


### Fit baseline data by a nonparametric local smoothing procedure
### (LOESS is used here) to obtain an estimator of f_24(x). To accommodate
### a possible jump at time=12 p.m., the model is fitted separately in
### time intervals [0,12] and (12,24].
	
bhbfit1 <- loess.smooth(time[1:144],bhbdat[1:144],span=0.3,
           family="gaussian",degree=2,evaluation=144)
bhbfit2 <- loess.smooth(time[145:288],bhbdat[145:288],span=0.3,
           family="gaussian",degree=2,evaluation=144)

### Plot baseline data and the fitted model.
	
par(mfg=c(1,1,2,2))
plot(time,bhbdat,xlab="Time (hours)",ylab="SLEEP",xaxt="n",yaxt="n",
     ylim=c(0,100),pch=".",mgp=c(2,1,0),cex=0.6)
lines(time[1:144],bhbfit1$y)
lines(time[145:288],bhbfit2$y)
points(time[1:144],rep(100,144),pch="-",cex=0.6)
points(time[145:288],rep(100,144),pch="o",cex=0.6)
axis(1,at=c(0,6,12,18,24),labels=c("0","6","12","18","24"),cex=0.6)
axis(2,at=c(0,25,50,75,100),labels=c("0","25","50","75","Light"),cex=0.6)
title(xlab="(a)",cex=0.7) 


### In order to fit the test data by a piecewise polynomial, we need to
### define a time variable which repeats time interval (0,6] four times,
### corresponding to four periods. An indicator variable "ind1" is also
### defined to accommodate a possible jump at the middle position of
### each period. 

time1 <- rep(time[1:72],4)
ind1  <- rep(c(rep(0,36),rep(1,36)),4)


### Fit the test data with the estimator of f_24(x) as a covariate
### variable. From this step (corresponding to model (2.1) in the paper),
### an estimator of f_6(x) is determined. Because polynomial regression
### is used in this step, the backward model selection procedure should
### be used with the hierarchy principle to determine the orders of the
### related polynomials.
	
bhtfit <- lm(bhtdat ~ c(bhbfit1$y,bhbfit2$y) + time1 + time1^2 + 
	     time1^3 + time1^4 + ind1 + ind1*time1 + ind1*time1^2 + 
	     ind1*time1^3 + ind1*time1^4)


### Plot the estimator of f_6(x).
	
par(mfg=c(2,1,2,2))
plot(time[1:72],bhtfit$fitted.values[1:72]-bhtfit$coef[2]*bhbfit1$y[1:72], 
     type="l",xlab="Time (hours)",ylab="SLEEP",xaxt="n",yaxt="n",
     ylim=c(0,100),xlim=c(0,24),mgp=c(2,1,0),cex=0.6)
lines(time[73:144],bhtfit$fitted.values[1:72]-bhtfit$coef[2]*bhbfit1$y[1:72])
lines(time[145:216],bhtfit$fitted.values[1:72]-bhtfit$coef[2]*bhbfit1$y[1:72])
lines(time[217:288],bhtfit$fitted.values[1:72]-bhtfit$coef[2]*bhbfit1$y[1:72])
points(time[1:36],rep(100,36),pch="-",cex=0.6)
points(time[37:72],rep(100,36),pch="o",cex=0.6)
points(time[73:108],rep(100,36),pch="-",cex=0.6)
points(time[109:144],rep(100,36),pch="o",cex=0.6)
points(time[145:180],rep(100,36),pch="-",cex=0.6)
points(time[181:216],rep(100,36),pch="o",cex=0.6)
points(time[217:252],rep(100,36),pch="-",cex=0.6)
points(time[253:288],rep(100,36),pch="o",cex=0.6)
axis(1,at=c(0,6,12,18,24),labels=c("0","6","12","18","24"),cex=0.6)
axis(2,at=c(0,25,50,75,100),labels=c("0","25","50","75","Light"),cex=0.6)
title(xlab="(c)",cex=0.7)


### Fit the residuals of model (2.1) by a nonparametric local smoothing
### procedure (LOESS is used here). From this step, an estimator of
### f_a(x) could be obtained (cf. model (2.2)).
	
resifit <- loess.smooth(time,bhtfit$residuals,span=0.3,
           family="gaussian",degree=2,evaluation=288)


### Plot the estimator of f_a(x).
	
par(mfg=c(2,2,2,2))
plot(time,bhtfit$residuals,xlab="Time (hours)",pch=".",ylab="SLEEP",
     xaxt="n",yaxt="n",ylim=c(-40,40),mgp=c(2,1,0),cex=0.6)
lines(time,resifit$y)
points(time[1:36],rep(40,36),pch="-",cex=0.6)
points(time[37:72],rep(40,36),pch="o",cex=0.6)
points(time[73:108],rep(40,36),pch="-",cex=0.6)
points(time[109:144],rep(40,36),pch="o",cex=0.6)
points(time[145:180],rep(40,36),pch="-",cex=0.6)
points(time[181:216],rep(40,36),pch="o",cex=0.6)
points(time[217:252],rep(40,36),pch="-",cex=0.6)
points(time[253:288],rep(40,36),pch="o",cex=0.6)
axis(1,at=c(0,6,12,18,24),labels=c("0","6","12","18","24"),cex=0.6)
axis(2,at=c(-40,-20,0,20,40),labels=c("-40","-20","0","20","Light"),cex=0.6)
title(xlab="(d)",cex=0.7)


### Combine the estimators of f_24(x), f_6(x) and f_a(x) to get a fitted
### model of (2.2). This fitted model and the test data are plotted as
### follows.
	
par(mfg=c(1,2,2,2))
plot(time,bhtdat,type="p",xlab="Time (hours)",pch=".",ylab="SLEEP",
     xaxt="n",yaxt="n",ylim=c(0,100),mgp=c(2,1,0),cex=0.6)
lines(time[1:72],bhtfit$fitted.values[1:72]+resifit$y[1:72])
lines(time[73:144],bhtfit$fitted.values[73:144]+resifit$y[73:144])
lines(time[145:216],bhtfit$fitted.values[145:216]+resifit$y[145:216])
lines(time[217:288],bhtfit$fitted.values[217:288]+resifit$y[217:288])
points(time[1:36],rep(100,36),pch="-",cex=0.6)
points(time[37:72],rep(100,36),pch="o",cex=0.6)
points(time[73:108],rep(100,36),pch="-",cex=0.6)
points(time[109:144],rep(100,36),pch="o",cex=0.6)
points(time[145:180],rep(100,36),pch="-",cex=0.6)
points(time[181:216],rep(100,36),pch="o",cex=0.6)
points(time[217:252],rep(100,36),pch="-",cex=0.6)
points(time[253:288],rep(100,36),pch="o",cex=0.6)
axis(1,at=c(0,6,12,18,24),labels=c("0","6","12","18","24"),cex=0.6)
axis(2,at=c(0,25,50,75,100),labels=c("0","25","50","75","Light"),cex=0.6)
title(xlab="(b)",cex=0.7)


graphics.off()