############################################################################
#						 		           #
#  Program to perform iteratively reweighted least squares by              #
#  by the simple Gauss-Newton approach (with no fancy modifications).      #
#								           #
#  The user must define 3 functions below:                                 #
#								           #
#	f(x,b)	the nonlinear regression model depending on a (n x k)      #
#		matrix of predictor variables and a (p x 1) regression     #
#		parameter b                                                #
#								           #
#	fb(x,b)	the (n x p) matrix of derivatives of f wrt b               #
#								           #
#	g(b,theta,z)  the variance model with variance parameters theta    #
#								           #
#  The user must also provide the following:                               #
#								           #
#	y	(n x 1) data vector                                        #
#								           #
#	x	(n x k) precictor values                                   #
#								           #
#	bstart	starting value for b                                       #
#								           #
#	theta	the fixed, known value for theta                           #
#								           #
#	tol	convergence criterion -- if 2 successive iterates of b     #
#		are within tol of each other relative to the current       #
#		value, declare that the algorithm has converged            #
#								           #
#	imax	maximum number of iterations to perform                    #
#								           #
############################################################################

#  logistic regression

#  Functions defining f, fb, g

#  Regression model is the logistic model

f <- function(x,b){
	temp <- exp(b[1] + b[2]*x[,1] + b[3]*x[,2] + b[4]*x[,3]
	 + b[5]*x[,4] + b[6]*x[,5] + b[7]*x[,6] + b[8]*x[,7])
	f1 <- temp/(1+temp)
	f1
}

fb <- function(x,b){
	n <- nrow(x)
	p <- length(b)
	fb1 <- matrix(0,n,p)
	fb1[,1] <- f(x,b)*(1-f(x,b))
	fb1[,2] <- f(x,b)*(1-f(x,b))*x[,1]
	fb1[,3] <- f(x,b)*(1-f(x,b))*x[,2]
	fb1[,4] <- f(x,b)*(1-f(x,b))*x[,3]
	fb1[,5] <- f(x,b)*(1-f(x,b))*x[,4]
	fb1[,6] <- f(x,b)*(1-f(x,b))*x[,5]
	fb1[,7] <- f(x,b)*(1-f(x,b))*x[,6]
	fb1[,8] <- f(x,b)*(1-f(x,b))*x[,7]
	fb1
}

g <- function(b,theta,z){ 

#  Bernoulli variance
	
	g1 <- sqrt(f(z,b)*(1-f(z,b)))
	g1
}

#  data

thedata <- scan("trial.dat")
thedata <- matrix(thedata,ncol=8,byrow=T)

#  need to create the 2 dummy variables for smoking

y <- thedata[,2]
covs <- thedata[,3:8]

x3 <- covs[,3]==1
x4 <- covs[,3]==2


n <- length(y)

# form the appropriate vector x as in problem

x <- cbind(covs[,1],covs[,2],x3,x4,covs[,4:6])

#  starting value for beta

xstar <- cbind(rep(1,n),x)
ystar <- (y+1/2)/2
wstar <- diag(ystar*(1-ystar))
eta <- log(ystar/(1-ystar))
zstar <- eta + (y-ystar)/(ystar*(1-ystar))

bstart <- solve(t(xstar)%*%wstar%*%xstar)%*%t(xstar)%*%wstar%*%zstar
	
#  convergence tolerance and max number of iterations
	
tol <- 1e-8
imax <- 500

#  name of output file
	
outfile <- "trial.irwls"

cat("FITS OF TRIAL DATA BY LOGISTIC REGRESSION",
	file=outfile,"\n","\n",append=F)

########################################################################

#  generic code to perform irwls and compute sigma

n <- length(y)
p <- length(bstart)

#  function specifying convergence criterion
	
converged <- function(tol,db,iter,imax){
	bmax <- max(abs(db))
	gmax <- iter==imax
	bmax <- bmax<tol
	converged <- gmax | bmax
	converged
}

#  function to perform IRWLS
	
irwls <- function(y,x,b,theta,tol,imax,f,fb,g){

#  initialize

	db <- 1
	iter <- 1
         

#  start iterating - stop when converged
	
	while(! converged(tol,db,iter,imax)){

#  current gradient matrix and weights
	
	  Xb <- fb(x,b)
	  w <- 1/g(b,theta,x)^2

#  construct X'WX
	
	  XbWXb <- t(Xb) %*% diag(w) %*% Xb

#  constructed vaariable z

	  z <- Xb %*% b + (y-f(x,b))
	  XbWz <- t(Xb) %*% diag(w) %*% z

#  get new iterate

          bnew <- solve(XbWXb) %*% XbWz
    
#  compute relative change


	  db <- (bnew-b)/b

	  b <- bnew
	  iter <- iter+1
  }


#  list of results returned by function
	
	list(beta=b,iend=iter-1)

}

theta <- 1  #  there really is no theta here

result <- irwls(y,x,bstart,theta,tol,imax,f,fb,g)
betahat <- result$beta

cat("(b) Starting value for beta = ",round(bstart,6),file=outfile,"\n","\n",
   append=T)

cat("(c) Estimate of beta = ",round(betahat,6),file=outfile,"\n",append=T)
cat("Number of iterations",result$iend,file=outfile,"\n","\n",append=T)

x60 <- c(0,0,0,0,60,85,130)
odds60 <- exp(betahat[1]+t(betahat[2:8])%*%x60)

cat("(e) Odds that a randomly chosen 60 year old who was randomized
to placebo and who had never smoked, with no previous MI, weighing 85
kg, and with systolic blood pressure of 130 will experience death or
MI in 30 days = ", odds60,file=outfile,"\n","\n",
   append=T)