#
# Epage > 137
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

if(is.null(version$language) == FALSE){ 
  require(alr3)
}else{
  library(alr3)
}

if(is.null(version$language) == FALSE) data(sleep1)

# local assignment required for use with Splus:
sleep <- sleep1
sleep$lB <- logb(sleep$BodyWt,2)
sleep$D <- factor(sleep$D)

# create the indicators for the categories of G for g = 1, 2, \cdots, 5 
sleep$D1 <- sleep$D2 <- sleep$D3 <- sleep$D4 <- sleep$D5 <- rep(0,dim(sleep)[1])
sleep$D1[sleep$D==1] <- 1
sleep$D2[sleep$D==2] <- 1
sleep$D3[sleep$D==3] <- 1
sleep$D4[sleep$D==4] <- 1
sleep$D5[sleep$D==5] <- 1

# to get the initial guess at this model ... take \gamma=0 and fit the
# regression model where we have different slopes for each of the various factors 
#
m1 <- lm(TS ~ lB:D, data=sleep, na.action=na.omit)

b0 <- coefficients(m1)[1]
s1 <- coefficients(m1)[1+1]
s2 <- coefficients(m1)[2+1]
s3 <- coefficients(m1)[3+1]
s4 <- coefficients(m1)[4+1]
s5 <- coefficients(m1)[5+1]

m2 <- nls(TS ~ beta0 + ( D1*beta11 + D2*beta12 + D3*beta13 + D4*beta14 + D5*beta15 ) * ( lB + gamma ),
          data=sleep, na.action=na.omit,
          start=list(beta0=b0,beta11=s1,beta12=s2,beta13=s3,beta14=s4,beta15=s5,gamma=2.5) )

attach(sleep)
plot( lB, TS )

# for plotting the different models on top of each other ... not tested or run 
if( 0 ) { 
  attach(tdata)
  plot(A,Gain,pch=S)
  xx <- seq(0,.44,length=100)
  # set m=1 in predict because of the weights
  lines(xx,predict(m2,data.frame(m=1,A=xx,S1=rep(1,100),
                                 S2=rep(0,100),S3=rep(0,100))),lty=1)
  lines(xx,predict(m2,data.frame(m=1,A=xx,S1=rep(0,100),
                                 S2=rep(1,100),S3=rep(0,100))),lty=2)
  lines(xx,predict(m2,data.frame(m=1,A=xx,S1=rep(0,100),
                                 S2=rep(0,100),S3=rep(1,100))),lty=3)
}