#############################################
# HW 7, Splus code for question 1
#############################################

# start graphics window
motif()
par(mfrow=c(2,2),cex=0.7,mex=0.7,mar=c(4,4,0.1,0.1),mgp=c(0.7,0,0))

read.dta <- function()
{
  dta <<- scan(skip=2,
	       file="214.dta",what=list(patid=0,y=0,ctx=0,gm=0,day=0))
  x <<- dta[[5]]
  ctx <<- dta[[3]]
  gm <<- dta[[4]]
  y <<- dta[[2]]
  patient <<- dta[[1]]
}

# loop over all patients, plot data for i-th patient
# and compute m.l.e. for logistic fit
# return matrix with logistic regr par:
#     (Nx2) matrix b
# and covariate matrices:  (NxT) X, (Nx2): X2=[ctx,gm]
#     matrix of responses: (NxT) Y
# where T=12 is the max number responses per patient
# missing values are filled with 0

plot.all <- function()
{ 
  b <<- NULL # will record m.l.e. for each patient in i-th row
  ni <<-NULL # will record number obs per patient
  day <<- matrix(0,nrow=52,ncol=12)
  X <<- matrix(0,nrow=52,ncol=2)
  Y <<- matrix(0,nrow=52,ncol=12)
  i _ 0
  repeat{
    i _ i + 1
    idx _ patient==i
    n <- sum(idx)
    ni <<- c(ni,n)
    day[i,1:n] <<- x[idx]
    X[i,1] <<- ctx[idx][1]
    X[i,2] <<- gm[idx][1]
    Y[i,1:n] <<- y[idx]
    if (sum(idx)<7)
      b <<- rbind(b,b[nrow(b),])
    else{
      plot(x[idx],y[idx],xlab="DAY",ylab="WBC")
      title(sub=paste("PATIENT",i))
      g _ glim(x[idx],y[idx],error="gaussian", iter.max=100,
	       link="logit",intercept=T, init = c(-0.7,0.55))
      ehat <- x[idx]*g$coef[2]+g$coef[1]
      yhat <- 1/(1+exp(-ehat))
      lines(x[idx],yhat)
      b <<- rbind(b,g$coef)
      cat("i=",i,"\t b=",format(b[i,]),"\t n=",sum(idx),"\n")
      if (i==52) break
    }
  }
  mb _ apply(b[-c(1,15,46,47),],2,mean)
  b[1,] <<- mb
  b[15,] <<- mb
  b[46,] <<- mb
  b[47,] <<- mb
  NULL # values are returned in global variables b,ni,day,X,Y
}


init.bugs <- function()
{ # prepares data files for bugs
  options(digits=3)
  write(format(t(X)),file="X",ncolumns=ncol(X))
  write(format(t(Y)),file="Y",ncolumns=ncol(Y))
  write(format(t(day)),file="day",ncolumns=ncol(day))
  write(format(ni),file="ni",ncolumns=1)
  write(format(t(b)),file="b",ncolumns=ncol(b))
}