REML = function(tol=1e-6,maxit=30,data=s$data,
                    model=4){

logHR = NULL ; se.logHR = NULL
M     = length(table(data$trial))
N     = dim(data)[1]
Z     = matrix(nrow=N,ncol=M); Z[,] = 0
Qsum  = 1 ; count = 0  ; code  = 1

for(i in 1:M){
time          = data$time[data$trial==i]
event         = data$event[data$trial==i]
x             = data$x[data$trial==i]
fit           = coxph(Surv(time,event) ~ x)
logHR[i]      = fit$coef
se.logHR[i]   = summary(fit)$coef[3]

}

# Starting values

y = logHR ; k = length(y); s = se.logHR ; w = 1/s^2 ; sum.w  = sum(w)
mu.hat  = sum(y*w)/sum.w  ; W = sum.w - (sum(w^2)/sum.w)
Q = sum(w*(y-mu.hat)^2) ; tausq.hat  = max(0,(Q-(k-1))/W) 
w.tauhat = 1/(s^2+tausq.hat); mu.tau.hat = sum(w.tauhat*y)/sum(w.tauhat)

u      = matrix(nrow=M,ncol=1)
Lambda = se.logHR^2/(se.logHR^2+tausq.hat)
u[,]   = Lambda*mu.tau.hat + (1-Lambda)*logHR
beta   = matrix(nrow=(M-1),ncol=1); beta[,] = 0
theta  = matrix(nrow=1,ncol=1)  ; theta[,]  = mu.tau.hat
tausq  = max(mu.tau.hat,0.05)

if(model==4){

data1 = data[order(data$trial,data$time),]

print("fitting stratified model:")
print("combined parameter change")

dim0  = table(data1$trial) ; dimU = cumsum(c(0,dim0))
delta = as.matrix(data1$event,ncol=1)
X     = as.matrix(data1$x,ncol=1)
for(jj in 1:N){Z[jj,data1$trial[jj]]=X[jj]} ; XZ = cbind(X,Z)
P = rbind(theta,u)

while(Qsum >= tol){

# sum first likelihood term
# over each trial

newU        = P[2:(M+1)]
count       = count + 1
dldeta      = matrix(nrow = N,ncol=1)
d2ldeta2    = matrix(nrow = N,ncol=N)
d2ldeta2[,] = 0
Tdiag       = diag(c(0,rep(1/tausq,M)))
Uvec        = as.matrix(c(0,newU/tausq))

for(jj in 1:M){

low      = dimU[jj]+1
up       = dimU[jj+1]
delta_st = as.matrix(delta[low:up],ncol=1)
X_st     = matrix(X[low:up,],ncol=1)
Z_st     = matrix(Z[low:up,jj],ncol=1)
XZ_st    = cbind(X_st,Z_st)
R_st     = matrix(nrow = dim0[jj],ncol=dim0[jj])
R_st[,]  = 1 ; R_st[lower.tri(R_st)]=0

TAB  = as.vector(table(data1$time[data1$trial==jj]))
CTAB = cumsum(TAB)
for(gg in 1:length(TAB)){
z     = TAB[gg]
if(z != 1){
LOW    = CTAB[gg] - TAB[gg] +1 ; UP  = CTAB[gg]
Dim    = UP-LOW+1 ;   Newmat = diag(Dim)
for(kk in 1:Dim){Newmat[kk,] = rep((Dim-kk+1)/Dim,Dim)}
R_st[LOW:UP,LOW:UP] = Newmat

}}


P_st       = rbind(P[1],newU[jj])
eta        = XZ_st%*%P_st
e.eta      = exp(eta)
SUMe.eta   = R_st%*%e.eta
W          = diag(e.eta[,])
a          = delta_st/SUMe.eta
A          = diag(a[,])
Mmat       = diag(rep(1,dim0[jj])); Mmat[lower.tri(Mmat)]=1
b          = cumsum(a)
B          = diag(b)
Ones       = as.matrix(rep(1,dim0[jj]),ncol=1)

dldeta[low:up,1]                = delta_st - W%*%Mmat%*%A%*%Ones
d2ldeta2[low:up,low:up] = W%*%B - W%*%Mmat%*%(A^2)%*%t(Mmat)%*%W

}

# calculate Pnew
# given tausq

V            = t(XZ)%*%d2ldeta2%*%XZ + Tdiag
Vinv         = solve(V)
Pnew         = P + Vinv%*%t(XZ)%*%dldeta - Vinv%*%Uvec
newU         = Pnew[2:(M+1)]
REMLtrace    = sum(diag(Vinv)[2:(M+1)])
tausqnew     = tausq*(t(newU)%*%newU)/(M*tausq -REMLtrace)

# update

Qsum  = abs(sum(P-Pnew))
P     = Pnew
tausq = tausqnew
print(Qsum)

if(count==30){Qsum = 0 ; print("max.iteration reached") ; code = 0 }

}
}


if(model==5){

print("fitting fixed trial random treatment effect model:")
print("combined parameter change")

data1 = data[order(data$time),]
dim0  = table(data1$trial)
X     = matrix(nrow=N,ncol=M); X[,] = 0 ; X[,1] = data1$x
R     = matrix(nrow=N,ncol=N)      ; R[,] = 1 ; R[lower.tri(R)]=0
tBH   = rep(1,N)                   ; tBH[data1$trial==1] = 0
delta = as.matrix(data1$event,ncol=1)

for(jj in 1:N){
temp           = data1$trial[jj]
Z[jj,temp]     = X[jj,1]
if(tBH[jj]==1){X[jj,temp] = tBH[jj]}
}

XZ  = cbind(X,Z) ; P = rbind(theta,beta,u)

###################
# Deal with Ties  #
###################

TAB  = as.vector(table(data1$time))
CTAB = cumsum(TAB)
for(gg in 1:length(TAB)){
z     = TAB[gg]
if(z != 1){
LOW    = CTAB[gg] - TAB[gg] +1 ; UP  = CTAB[gg]
Dim    = UP-LOW+1 ;   Newmat = diag(Dim)
for(jj in 1:Dim){Newmat[jj,] = rep((Dim-jj+1)/Dim,Dim)}
R[LOW:UP,LOW:UP] = Newmat
}}

while(Qsum >= tol){

count = count + 1 ; eta  = XZ%*%P
e.eta = exp(eta)  ; SUMe.eta  = R%*%e.eta
W     = diag(e.eta[,]); a = delta/SUMe.eta
A     = diag(a[,]);  b  = cumsum(a) ; B = diag(b)
Mmat  = diag(rep(1,N)); Mmat[lower.tri(Mmat)]=1
Ones  = as.matrix(rep(1,N),ncol=1)

newU       = P[(M+1):(2*M)]
Tdiag      = diag(c(rep(0,M),rep(1/tausq,M)))
Uvec       = as.matrix(c(rep(0,M),newU/tausq))
dldeta     = delta - W%*%Mmat%*%A%*%Ones
d2ldeta2   = W%*%B - W%*%Mmat%*%(A^2)%*%t(Mmat)%*%W
V          = t(XZ)%*%d2ldeta2%*%XZ + Tdiag
Vinv       = solve(V)
Pnew       = P + Vinv%*%t(XZ)%*%dldeta - Vinv%*%Uvec
newU       = Pnew[(M+1):(2*M)]
REMLtrace  = sum(diag(Vinv)[(M+1):(2*M)])
tausqnew   = tausq*(t(newU)%*%newU)/(M*tausq -REMLtrace)
Qsum       = abs(sum(P-Pnew))
print(Qsum)
P          = Pnew ; tausq = tausqnew

if(count==30){Qsum = 0 ; print("max.iteration reached") ; code = 0 }

}}

var.theta = Vinv[1,1]
CI.theta  = P[1,1] + c(-1,1)*1.96*sqrt(var.theta)

return(list(tausq=tausq,P=P,count=count,code=code,
            var.theta=var.theta,CI.theta=CI.theta))

}