R/lognorm1.iter.R

In frailtyHL: Frailty Models via Hierarchical Likelihood

lognorm1.iter <- function(X,B,Z,v,penalty,alpha,tun1=NULL,tun2=NULL,di,Mi,idx2,del,varfixed){

n<-nrow(X)
p<-ncol(X)


iB<-B
iv<-v
ialpha<-alpha

rep<-1
err<-3
convergeind<-c(0)

while((rep<2000)&(err>0.0000001)&(err<5)){

if (varfixed==TRUE) result<-lognorm(X=X,B=iB,Z=Z,v=iv,penalty=penalty,alpha=ialpha,tun1=tun1,tun2=tun2,di,Mi,idx2,del,varfixed)
else result<-lognorm1(X=X,B=iB,Z=Z,v=iv,penalty=penalty,alpha=ialpha,tun1=tun1,tun2=tun2,di,Mi,idx2,del,varfixed)

err<-sum(abs(ialpha-result[[3]]))

if (varfixed==TRUE){err<-max(abs(iB-result[[1]]))}

#print("==== rep & err ===")
#print(c(rep,err))
#print(ialpha)

if (err<=0.0000001) {convergeind<-c(1)} else {convergeind<-c(-1)}

iB<-result[[1]]
iv<-result[[2]]
ialpha<-result[[3]]
mat<-result[[4]]
WL<-result[[5]]
Bi<-result[[6]]
hft<-result[[7]]
rep<-rep+1

#print(ialpha)
#print(convergeind)
#print(iB)
#print(iv)
}

#print(iB)

Zs<-Z
nrand<-length(Zs)
qs<-rep(0,nrand)

z<-NULL;v_h1<-NULL; alpha_h1<-NULL; rand.idx<-NULL

for (i in 1:nrand){
        z<-cbind(z,Zs[[i]])
        qs[i]<-dim(Zs[[i]])[2]
        v_h1<-c(v_h1,iv[[i]])
        alpha_h1<-c(alpha_h1,rep(ialpha[i],qs[i]))
        rand.idx<-c(rand.idx, rep(i, qs[i]))    
}

q<-sum(qs)

beta_h1<-iB

x<-X

mat<-mat

u_h1=exp(v_h1)

mat11=t(x)%*%mat%*%x
mat12=t(x)%*%mat%*%z
mat13=t(z)%*%mat%*%z

U=(1/alpha_h1)*diag(rep(1,q))
#print(det(mat13+U))

if ((err>5)|(det(mat13+U)<0.0000001)){BIC<-100000; se_beSp<-rep(0,len=p)}else{mmat=mat11-mat12%*%solve(mat13+U)%*%t(mat12)
hminv=solve(mmat+n*WL)
se_bem=sqrt(diag(hminv))
t=beta_h1/se_bem

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

eta=as.numeric(x%*%beta_h1 + z%*%v_h1)
expeta=exp(eta)
Wi=diag(expeta)

minv=solve(-mmat-n*WL)
cov=mmat
iHs=(-minv)%*%cov%*%(-minv)
se_beS=sqrt(diag(iHs))
se_beSp=se_beS[1:p]
t1=beta_h1/se_beSp

eta=as.numeric(x%*%beta_h1 + z%*%v_h1)
expeta=exp(eta)
one=rep(1,len=n)
done=rep(1,len=idx2)
oq=rep(1,q)

pi=3.14159265359
term0=as.numeric(t(Mi)%*%expeta)

hlike1= (t(one)%*%(del*eta))-(t(done)%*%(di*log(term0)))
hlike2=-(q/2)*log(2*pi)-(1/2)*sum(log(alpha_h1))
hlike3=-0.5*(t(v_h1)%*% diag(1/alpha_h1) %*%v_h1)

hlike=hlike1+hlike2+hlike3

mat0=(Wi%*%Bi)
H22=t(z)%*%mat%*%z+U

oq=rep(1,len=q)
eiv=eigen(H22)$value
adj2=as.numeric(-0.5*(t(oq)%*%log(eiv/(2*pi))))
pvh=hlike+adj2

if(length(beta_h1)!=1){
WL=diag(hft/(abs(beta_h1)))
}else{
WL=hft/(abs(beta_h1))
}
iH=solve(mmat+n*WL)%*%(mmat)
elam=sum(diag(iH))

########################################
if (varfixed==TRUE) {pvh<-hlike1}
########################################
 
 BIC = -2*pvh+log(n)*elam
}

if((penalty=="lasso")|(penalty=="scad")){
result<-list(round(iB,digits=5),round(se_beSp,digits=5),round(v_h1,digits=5),round(ialpha,digits=5),tun1,round(BIC,digits=5),rep,convergeind)
names(result)<-c("beta","se","v","alpha","tuning_parameter","BIC","repetition","convergence")
return(result)
}
if(penalty=="hl"){
result<-list(round(iB,digits=5),round(se_beSp,digits=5),round(v_h1,digits=5),round(ialpha,digits=5),c(tun1,tun2),round(BIC,digits=5),rep,convergeind)
names(result)<-c("beta","se","v","alpha","tuning_parameter","BIC","repetition","convergence")
return(result)
}

}