#  SPLUS PROGRAM FOR MAXIMUM LIKELIHOOD (ML) LINEAR REGRESSION IN 
#  THE PRESENCE OF TYPE II LEFT AND INTERVAL CENSORING OF THE 
#  OUTCOME VARIABLE. IT IS ASSUMED THAT SOME OBSERVATIONS ARE 
#  "NON-DETECTS", I.E. THEY FALL BELOW THE DETECTION LEVEL OF THE
#  INSTRUMENT AND ARE LEFT CENSORED, OTHERS ARE "TRACE" 
#  MEASUREMENTS, I.E. THEY LIE ABOVE THE LIMIT OF DETECTION, BUT 
#  ARE BELOW SOME LEVEL WHERE RELIABLE QUANTIFICATION IS 
#  POSSIBLE.  THESE OBSERVATIONS ARE INTERVAL CENSORED.
#  ALTHOUGH THIS PROBLEM IS FRAMED IN TERMS OF LEFT CENSORING, 
#  THE METHODOLOGY APPLIES EQUALLY TO THE RIGHT CENSORING 
#  SETTING, SIMPLY BY LOOKING AT THE MIRROR IMAGE.

# 29 November, 1999

#  THIS IS A MAXIMUM-LIKELIHOOD TECHNIQUE THAT WILL ESTIMATE THE
#  COEFFICIENTS AND VARIANCE USING AN EM ALGORITHM
#  TO USE THIS FUNCTION, THE FOLLOWING ARE REQUIRED:
#	- VECTOR OF OUTCOMES, y; (n x 1)
#	- IT IS ASSUMED THAT y IS ALREADY TRANSFORMED TO BE
#  	  APPROXIMATELY NORMALLY DISTRIBUTED.
#	- CENSORED VALUES FOR y ARE REPLACED BY SENSIBLE FIRST
#       ESTIMATES
#       FOR EXAMPLE, THE MIDPOINT OF THE RELEVANT CENSORING INTERVAL
#	- MATRIX OF DEPENDENT VARIABLES, x, NOT INCLUDING
#       INTERCEPT;(n x p)
#	- VECTOR OF "FLAGS" (flag) DESCRIBING FOR EACH OBSERVATION
#  	  WHETHER IT IS AT THE NON-DETECT ("ND")
#	  OR WHETHER IT IS AT TRACE ("<QL")
#	  OR NOT CENSORED AT ALL ("[no"); (n x 1)
#  A VECTOR CONTAINING THE NAMES OF THE X-VARIABLES IS ALSO TO BE
#  INCLUDED WHEN THE FUNCTION IS CALLED, FOR EXAMPLE,
#  names.vector<- c("age","weight") 
#  NOTE THAT LENGTH OF EACH NAME SHOULD BE LESS THAN 9 LETTERS
#  EACH ALSO NEED TO INCLUDE THE NAME OF THE OUTCOME (DEPENDENT)
#  VARIABLE FOR EXAMPLE, y.name<- "DMTP"



maxlikelihd.function<-function(x, y, flag, ND.limit, QL.limit, difference=0.01, xnames.vector,y.name)
{

#  AN INTERNAL FUNCTION FOR PRINTING AT A LATER STAGE
my.format <- function (S, len, justify="l", ...)
{
        S2 <- format(S,...)
        x <- "                             "
        n <- len-nchar(S2)
        if (justify=="l")
                paste(S2,substring(x,1,n),sep="")
        else
                paste(substring(x,1,n),S2,sep="")
}

sample.l	<-length(y)
marker	<-0

this.dataframe<-data.frame(x=x, y=y, flag=flag, ND.limit=ND.limit, QL.limit=QL.limit)
attach(this.dataframe)

# PROGRAM:
# THIS PROGRAM FITS A LINEAR REGRESSION MODEL TO LEFT-CENSORED 
# INTERVAL-CENSORED DATA VIA MAXIMUM LIKELIHOOD ESTIMATION (MLE)

temp	<- y

#  TO GET INITIAL ESTIMATES OF THE BETAS (old.betas)
#  BASED ON INITIAL ESTIMATE OF y.

lmmodel1		<- lm(y ~ x)
lmmodel1.sum	<- summary(lmmodel1)
old.betas		<- (lmmodel1$coefficients)

#  OBTAIN MODEL PARAMETERS BASED ON INITIAL ESTIMATE OF Y

lmmodel.temp1		<- lm(temp ~ x)
lmmodel.temp1sum		<- summary(lmmodel.temp1)
mlbetas.initial		<- lmmodel1$coefficients
mlvariance.initial	<- (lmmodel1.sum$sigma)^2
tvalues.initial		<- lmmodel1.sum$coefficients[2,3] 

sigmasq		<-(lmmodel1.sum$sigma)^2
sqrtsigmasq		<-sqrt(sigmasq)

#  SET UP THE X-MATRIX 
#  n IS THE NUMBER OF OBSERVATIONS IN THE DATASET.

n		<-length(y)  
ones.vector	<-as.vector(rep(1,n))
x.variables	<-as.matrix(cbind(ones.vector,x))

# INITIALISE MEAN
mew		<-as.vector(rep(0,n))
mew 		<-(x.variables %*% old.betas)

#  REPEATED ITERATIONS OF EM, WITH FILLING IN OF CENSORED Y
#  OBSERVATIONS BASED ON PREVIOUS ITERATIONS OF FITTED MODEL.
#  NOTE THAT CONVERGENCE IS DEFINED TO BE A MAXIMUM CHANGE IN THE
#  COEFFICIENTS OF LESS THAN 1%

yi.new		<-as.vector(rep(0,n))

numiterations	<-100

i			<-1
convergence		<-"F"

# LOOP
while ((convergence=="F") & (i<numiterations))
{

# CALCULATIONS TOWARDS DENOMINATOR OF MLE VARIANCE ESTIMATE
u<-as.vector((ND.limit-mew)/sqrt(sigmasq))
v<-as.vector((QL.limit-mew)/sqrt(sigmasq))

flag1.NDL	<-as.vector(ifelse(flag=="ND",1,0))

Dpart1<-as.vector(rep(0,n))
dnormu<-as.vector(dnorm(u))
pnormu<-as.vector(pnorm(u))

Dpart2<-as.vector(rep(0,n))
dnormv<-as.vector(dnorm(v))
pnormv<-as.vector(pnorm(v))

flag1.QL<-as.vector(ifelse(flag=="<QL",1,0))

for (i in 1:n)
{ Dpart1[i]<-(flag1.NDL[i]*dnormu[i]*u[i])/pnormu[i] 
  Dpart2[i]<-(flag1.QL[i]*(dnormv[i]*v[i]-dnormu[i]*u[i])/
  (pnormv[i]-pnormu[i])) 
 }

Dpart1.sum<-sum(Dpart1,na.rm=T)
Dpart2.sum<-sum(Dpart2,na.rm=T)

flag1.uncens<-as.vector(ifelse(flag=="[no",1,0))
num.uncensored<-sum(flag1.uncens)
D<-num.uncensored + Dpart1.sum + Dpart2.sum

uncens <-flag1.uncens*y

#  UPDATE NEW Y-VALUES FOR THE ND and <QL FLAGGED VALUES
#  NOTE THAT THE UNCENSORED Y-VALUES ARE NOT CHANGED

  for (i in 1:n)
  {
   if (flag[i]=="ND")  
   {yi.new[i]<-(mew[i]-(sqrtsigmasq*(dnormu[i]/pnormu[i]))) }

   if (flag[i]=="<QL")
   {yi.new[i]<-  (mew[i]-sqrtsigmasq*  ((dnormv[i]-dnormu[i])/
    (pnormv[i]-pnormu[i]))) }

   if (flag[i]=="[no")
   {yi.new[i]<-y[i]}
  }

#  FIT A LINEAR REGRESSION MODEL USING THE UPDATED DATA
lmmodel1		<- lm(yi.new~ x)
lmmodel1.sum	<-summary(lmmodel1)
new.betas		<-(lmmodel1$coefficients)

#  COMPARE THE COEFFICIENTS WITH THOSE FROM THE PREVIOUS ITERATION.
if (max((new.betas - old.betas)/old.betas)< difference)
{convergence<-"T"}

if (convergence=="F")
  {
     i<-i+1
     old.betas<-new.betas
  }

# RECALCULATE MEAN IF NOT YET CONVERGED
mew <-(x.variables %*% old.betas)

     #  CALCULATION OF ESTIMATE OF RESIDUAL VARIANCE
     
     mew.uncens<-as.vector(rep(0,n))

     for (j in 1:n)
       {
        if (abs(uncens[j])>0)  mew.uncens[j]<-mew[j]
        if (uncens[j]==0)      mew.uncens[j]<-0
       }

#  THE FINAL CALCULATION OF SIGMA-SQUARED     
     sigmasq			<-sum((uncens-mew.uncens)^2)/D
     sqrtsigmasq			<-sqrt(sigmasq)
     mlvariance.keep		<-sigmasq
   

}   #  END OF LOOP

mlbetas.keep	<-lmmodel1$coefficients
mlvariance.keep	<-sigmasq

#  CALCULATING THE INFORMATION MATRIX

sigma	<-sqrtsigmasq
c1	<-as.vector(ND.limit - mew)
c2	<-as.vector(QL.limit - mew)
phi1	<-dnorm(c1/sigma)
cphi1	<-pnorm(c1/sigma)
phi2	<-dnorm(c2/sigma)
cphi2	<-pnorm(c2/sigma)

numxvars<-dim(x.variables)[2]
inf.matrix<-matrix(0,ncol=(numxvars+1), nrow=(numxvars+1))

for(i in 1:n)
{

if (flag[i]=="ND")
{
#  CODE FOR THE ELEMENTS OF THE INFORMATION MATRIX TOP LEFT
#  POSITIONS
for (j in 1:numxvars)
  {
    for (k in j:numxvars)
    {
   inf.matrix[j,k]<-inf.matrix[j,k]+(c1[i]*phi1[i]*cphi1[i]/sigmasq+phi1[i]*
 phi1[i]/sigma)*x.variables[i,k]*x.variables[i,j]/cphi1[i]/cphi1[i]
    }
  }

#  CODE FOR THE BOTTOM ELEMENTS OF THE INFORMATION MATRIX
for (j in 1:numxvars)
  {
   inf.matrix[(numxvars+1),j]<- inf.matrix[(numxvars+1),j]+phi1[i]*(1- c1[i]^2
   /sigmasq-phi1[i]*c1[i]/sigma/cphi1[i])*x.variables[i,j]/cphi1[i]
 }

#  CODE FOR THE SINGLE RHS DIAGONAL BOTTOM ELEMENT
inf.matrix[numxvars+1,numxvars+1]<- inf.matrix[numxvars+1,numxvars+1] + phi1[i]*c1[i]*(3-(c1[i])^2/sigmasq-phi1[i]*c1[i]/sigma/cphi1[i])/4/cphi1[i]

} # END OF CODE FOR "ND" OBSERVATIONS

 
if (flag[i]=="<QL")
{
#  CODE FOR THE ELEMENTS OF THE INFORMATION MATRIX TOP LEFT
#  POSITIONS

for (j in 1:numxvars)
  {
    for (k in j:numxvars)
    {
    inf.matrix[j,k]<-inf.matrix[j,k]+(c2[i]*phi2[i]/sigmasq-
    c1[i]*phi1[i]/sigmasq)* x.variables[i,k]*x.variables[i,j]/(cphi2[i]-
    cphi1[i])+ x.variables[i,k]*x.variables[i,j]*(phi2[i]-phi1[i])^2 /(cphi2[i]-
    cphi1[i])^2/sigma
    }
  }

#  CODE FOR THE BOTTOM ELEMENTS OF THE INFORMATION MATRIX AND
#  CODE FOR THE RHS ELEMENTS OF THE INFORMATION MATRIX 
#  RESPECTIVELY

for (j in 1:numxvars)
  {
   inf.matrix[(numxvars+1),j]<- inf.matrix[(numxvars+1),j]+ (phi2[i]- phi1[i])*x.variables[i,j]/(cphi2[i]-cphi1[i])
   - x.variables[i,j]*(phi2[i]*c2[i]^2-phi1[i]*c1[i]^2)  /(cphi2[i]-cphi1[i])/sigmasq-(phi2[i]-phi1[i])*(phi2[i]*c2[i]-phi1[i]*c1[i])* x.variables[i,j]/(cphi2[i]-cphi1[i])^2/sigma

  }

#  CODE FOR THE SINGLE RHS DIAGONAL BOTTOM ELEMENT
inf.matrix[numxvars+1,numxvars+1]<- inf.matrix[numxvars+1,numxvars+1] +(phi2[i]*c2[i]-phi1[i]*c1[i])/(cphi2[i]-cphi1[i])*0.75
-(phi2[i]*c2[i]^3-phi1[i]*c1[i]^3)/(cphi2[i]-cphi1[i])/4/sigmasq -
(phi2[i]*c2[i]-phi1[i]*c1[i])^2/(cphi2[i]-cphi1[i])^2/4/sigma

} # END OF CODE FOR THE "QL" OBSERVATIONS

 
if (flag[i]=="[no")
{

# CODE FOR THE LHS ENTRIES IN THE INFORMATION MATRIX
  for (j in 1:numxvars)
  {
    for (k in j:numxvars)
    {
    	inf.matrix[j,k]<-inf.matrix[j,k]+(1/sigma)* x.variables[i,j]*  	x.variables[i,k]
    }
  } 

#  CODE FOR THE BOTTOM ELEMENTS OF THE INFORMATION MATRIX
for (j in 1:numxvars)
  {
   inf.matrix[(numxvars+1),j]<- inf.matrix[(numxvars+1),j]-2*(yi.new[i]-
   mew[i])*x.variables[i,j]/sigma  }

# CODE FOR THE SINGLE RHS DIAGONAL BOTTOM ELEMENT
inf.matrix[numxvars+1,numxvars+1]<- inf.matrix[numxvars+1,numxvars+1] +sigma/2-(yi.new[i]-mew[i])^2/sigma 

 } # END OF CODE FOR "[no" OBSERVATIONS

} # END OF CODE FOR ALL OBSERVATIONS

# TO "FILL IN" THE SYMMETRIC MEMBERS OF THE INFORMATION MATRIX
for (j in 1:numxvars)
  {
    for (k in j:numxvars)
    {
      inf.matrix[k,j]<-inf.matrix[j,k]
    }
  } 


#  FINAL TIDYING UP OF THE INFORMATION MATRIX 
#  CODE FOR THE TOP LHS OF MATRIX
for (i in 1:numxvars)
  {
    for (j in 1:numxvars)
    {
    inf.matrix[i,j]<- -inf.matrix[i,j]/sigma
    }
  }

#  CODE FOR THE BOTTOM ELEMENTS AND RHS ELEMENTS OF THE 
#  INFORMATION MATRIX 

for (i in 1:numxvars)
  {
  inf.matrix[(numxvars+1),i]<-inf.matrix[(numxvars+1),i]/sigma/sigmasq/2
   inf.matrix[i,(numxvars+1)]<- inf.matrix[i,(numxvars+1)]/sigma/sigmasq/2
  }

#  BOTTOM FINAL DIAGONAL ELEMENT OF THE INFORMATION MATRIX
inf.matrix[numxvars+1,numxvars+1]<-inf.matrix[numxvars+1,numxvars+1]/sigma/sigmasq/sigmasq
 
info.matrix		<- -inf.matrix
info.matrix		<-solve(info.matrix)

beta.variance	<-as.vector(rep(0,numxvars))

for (i in 1:numxvars)
{
beta.variance[i]	<-info.matrix[i,i]
}

#  NEW CODE TO ADJUST FOR A NEGATIVE BETA.VARIANCE
#  OVER EACH OF THE X-VARIABLES SEPARATELY

marker		<-as.vector(rep(0,numxvars))
beta.variancem	<-matrix(0,ncol=numxvars,nrow=numxvars)

for (i in 1:numxvars)
{
   if (beta.variance[i] <=0)
   {
     marker[i]<-1

     # Over all the observations in the dataset
       beta.variancem	<- sigmasq*solve(t(x.variables[flag=="[no",]) %*% 
 	 x.variables[flag=="[no",])
       beta.variance[i]	<-beta.variancem[i,i]
   }
}


tvalues.final <-new.betas/(sqrt(beta.variance))

percent.NDcen<-sum(ifelse(flag=="ND",1,0))/length(y)*100
percent.QLcen<-sum(ifelse(flag=="<QL",1,0))/length(y)*100
percent.totcen<-percent.NDcen + percent.QLcen


print("CENSORING DETAILS:")
print("------------------")
cat(paste("\n"))
cat(paste("Percentage ND censored:        ",round(percent.NDcen,1), "%","\n"))
cat(paste("Percentage <QL censored:       ",round(percent.QLcen,1), "%","\n"))
cat(paste("Total Percentage censored:     ",round(percent.totcen,1), "%","\n"))
cat(paste("\n"))

print("TABLE OF RESULTS - Maximum Likelihood")
cat(paste("\n"))
cat(paste("Outcome variable: ",y.name, "\n"))
cat(paste("\n"))
cat(paste("          Coefficient"," Std error ", "z-value"," P-value","\n"))
cat(paste("          -----------","","---------","", "-------","","--------","\n"))
cat(paste("(Initial)"))
cat(paste("\n"))
cat(
	my.format("Intercept",len=10, justify="l"),
	my.format(round(mlbetas.initial[1],4),len=9,justify="r"),
	my.format(round(lmmodel.temp1sum$coefficients[1,2],4),len=9	,justify="r"),
	my.format(round(lmmodel.temp1sum$coefficients[1,3],4),len=9	,justify="r"),
	my.format(round(lmmodel.temp1sum$coefficients[1,4],4),len=9	,justify="r"),
      "\n"
	)

for (i in 2:numxvars)
{
  cat(
	my.format(xnames.vector[i-1],len=10, justify="l"),
	my.format(round(mlbetas.initial[i],4),len=9,justify="r"),
	my.format(round(lmmodel.temp1sum$coefficients[i,2],4),len=9	,justify="r"),
	my.format(round(lmmodel.temp1sum$coefficients[i,3],4),len=9	,justify="r"),
	my.format(round(lmmodel.temp1sum$coefficients[i,4],4),len=9	,justify="r"),
	"\n"
     )
}

cat(paste("\n"))
cat(paste("(Final)"))
cat(paste("\n"))

cat(	
	my.format("Intercept",len=10, justify="l"),
	my.format(round(mlbetas.keep[1],4),len=9,justify="r"),
	my.format(round(sqrt(beta.variance[1]),4),len=9,justify="r"	),
	my.format(round(tvalues.final[1],4),len=9,justify="r"),
	my.format(round(2*(1-	pnorm(abs(tvalues.final[1]))),4),len=9,justify="r"),
      "\n"
	)

for (i in 2:numxvars)
{
	cat(
	my.format(xnames.vector[i-1],len=10, justify="l"),
	my.format(round(mlbetas.keep[i],4),len=9,justify="r"),
	my.format(round(sqrt(beta.variance[i]),4),len=9,justify="r"	),
	my.format(round(tvalues.final[i],4),len=9,justify="r"),
	my.format(round(2*(1-	pnorm(abs(tvalues.final[i]))),4),len=9,justify="r"),
	"\n"
      )
}

cat(paste("\n"))
cat(	my.format("Initial Residual Variance:",len=30, 	justify="l"),
	my.format(round(mlvariance.initial[1],4),len=8,justify="r	")
   )

cat(paste("\n"))
    
cat(	my.format("Final Residual Variance:",len=30, justify="l"),
	my.format(round(mlvariance.keep[1],4),len=8,justify="r")
   )
cat(paste("\n"))

if (sum(marker)>=1) {cat(paste("* Note:  Information matrix had negative variance (due to censoring) which has been replaced by a conservative estimate","\n"))} else if (sum(marker==0)) {cat(paste("\n"))}


#  A LIST OF POSSIBLE ITEMS THAT USERS MAY WISH TO ACCESS FROM
#  THE FUNCTION INCLUDING THE INITIAL AND FINAL BETA-COEFFICIENT
#  ESTIMATES, INITIAL AND FINAL VARIANCE ESTIMATES, R-SQUARED,
#  THE FITTED Y'S AND THE FINAL ESTIMATES OF THE Y-VALUES.	

list(	y.new				= yi.new, 		
	y.fitted			= lmmodel1$fitted.values, 
	mlbetas.initial		= mlbetas.initial, 
	mlbetas.final		= mlbetas.keep,
	R.squared			= lmmodel1.sum$r.squared, 
	mlvariance.initial	= mlvariance.initial, 
	mlvariance.final		= mlvariance.keep)

}  # END OF MAXIMUM LIKELIHOOD FUNCTION