Nothing
R/HGNlogL.R
In MixtureRegLTIC: Mixture Regression Models for Left-Truncated and Interval-Censored Data
Defines functions
LLF
GLLF
HLLF
HGNlogL
HGNlogLC
CMultiply
###################
### function: LLF()
###################
LLF=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype){
LLF=HGNlogL(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="LLF")
return(LLF)
}
####################
### function: GLLF()
####################
GLLF=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype){
GLLF=HGNlogL(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="GLLF")
return(GLLF)
}
####################
### function: HLLF()
####################
HLLF=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype){
HLLF=HGNlogL(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="HLLF")
return(HLLF)
}
#######################
### function: HGNlogL()
#######################
HGNlogL=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="LLF"){
HGNlogL=0
if (length(survtime)==6){
# Calculate truncation contribution
HGNlogL=HGNlogL-HGNlogLC(par.nfix,par.all,is.fix.par,survtime[4:6],x,is.reg,weight,mtype,outtype)
}
# Calculate interval censorship contribution
HGNlogL=HGNlogL+HGNlogLC(par.nfix,par.all,is.fix.par,survtime[1:3],x,is.reg,weight,mtype,outtype)
return(HGNlogL)
}
########################
### function: HGNlogLC()
########################
HGNlogLC=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="LLF"){
##### "par.all" => "par"
par=par.all
##### "par.nfix" => "par"
start=1
if (is.reg$beta){
len=length(par$beta); par$beta=par.nfix[seq(start,start+len-1)]; start=start+len
}
len=length(par$gamma); par$gamma=par.nfix[seq(start,start+len-1)]; start=start+len
len=length(par$alpha); par$alpha=par.nfix[seq(start,start+len-1)]; start=start+len
if (!is.fix.par$q){
len=length(par$q); par$q=par.nfix[seq(start,start+len-1)]; start=start+len
}
##### "y" and "nobs" (Survival Time and Number of Subjects)
y=survtime;nobs=dim(y)[1]
##### Mixture Type (The Type of Distribution Form for Each Subject)
ip2=rep(1,nobs) # "two components with cure"
index=which(mtype==1) # "one component"
if (length(index)>0){ip2[index]=0}
ip1=ip2 # used for derivative with beta
#######################################
##### "alpha" and "q" (scale and shape)
#######################################
alpha=A.fun(x$alpha,par$alpha)
q=Q.fun(x$q,par$q)
##########################################
##### Functions of "p1" & "p2=1-p1"
##### and It's First and Second Derivative
##########################################
##### "p1","p2" (function)
p1=rep(1,nobs)
if(is.reg$beta){
p1=p1.fun(x,par,is.reg)
index=which(mtype==1|mtype==3);if(length(index)>0) p1[index]=1
}
p2=rep(0,nobs);p2=1-p1
if(outtype=="GLLF"|outtype=="HLLF"){
##### "gp1","gp2" (first derivative)
gp1=list(beta=rep(0,nobs))
if (is.reg$beta){
gp1$beta=p1b.fun(x,par,is.reg)
index=which(mtype==1|mtype==3);if(length(index)>0) gp1$beta[index]=0
}
gp2=list(beta=rep(0,nobs));gp2$beta=-gp1$beta
if(outtype=="HLLF"){
##### "hp1", "hp2" (second derivative)
hp1=list(beta=rep(0,nobs))
if (is.reg$beta){
hp1$beta=p1bb.fun(x,par,is.reg)
index=which(mtype==1|mtype==3);if (length(index)>0) hp1$beta[index]=0
}
hp2=list(beta=rep(0,nobs));hp2$beta=-hp1$beta
}
}
################
##### Value of w
################
##### "w1","w2"
w1=rep(NA,nobs);w2=rep(NA,nobs)
index=which(y[,3]==1|y[,3]==3|y[,3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){w1[index]=w.fun(y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)}
index=which(y[,3]==3|y[,3]==4) # Interval or Left Censored Data
if (length(index)>0){w2[index]=w.fun(y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)}
##########################################
##### Density and Survival Functions of w
##### and It's First and Second Derivative
##########################################
##### Density and Survival Functions of w
### "fw1","Sw1","fw2","Sw2"
fw1=rep(NA,nobs);Sw1=rep(NA,nobs)
fw2=rep(NA,nobs);Sw2=rep(NA,nobs)
index=which(y[3]==1|y[3]==3|y[3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){
fw1[index]=dw.fun(w1[index],q[index])
Sw1[index]=Sw.fun(w1[index],q[index])
}
index=which(y[3]==3|y[3]==4) ### Interval Censored or Left Censored Data
if (length(index)>0){
fw2[index]=dw.fun(w2[index],q[index])
Sw2[index]=Sw.fun(w2[index],q[index])
}
if(outtype=="GLLF"|outtype=="HLLF"){
##### First Derivative of Density and Survival Functions of w
### "gfw1","gSw1","gfw2","gSw2"
gfw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
gSw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
gfw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
gSw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
index=which(y[3]==1|y[3]==3|y[3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){
gfw1$gamma[index]=dwg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gfw1$alpha[index]=dwa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){
gfw1$q[index]=dwq.fun(w1[index],q[index])
}
gSw1$gamma[index]=Swg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gSw1$alpha[index]=Swa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){gSw1$q[index]=Swq.fun(w1[index],q[index])}
}
index=which(y[3]==3|y[3]==4) # Interval Censored or Left Censored Data
if (length(index)>0){
gfw2$gamma[index]=dwg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gfw2$alpha[index]=dwa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){gfw2$q[index]=dwq.fun(w2[index],q[index])}
gSw2$gamma[index]=Swg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gSw2$alpha[index]=Swa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){gSw2$q[index]=Swq.fun(w2[index],q[index])}
}
if(outtype=="HLLF"){
##### Second Derivative of Density and Survival Functions of w
### "hfw1","hSw1","hfw2","hSw2"
hfw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
hSw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
hfw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
hSw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
index=which(y[3]==1|y[3]==3|y[3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){
hfw1$gamma[index]=dwgg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw1$ga[index]=dwga.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw1$gq[index]=dwgq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw1$alpha[index]=dwaa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw1$aq[index]=dwaq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw1$q[index]=dwqq.fun(w1[index],q[index])
hSw1$gamma[index]=Swgg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw1$ga[index]=Swga.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw1$gq[index]=Swgq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw1$alpha[index]=Swaa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw1$aq[index]=Swaq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw1$q[index]=Swqq.fun(w1[index],q[index])
}
index=which(y[3]==3|y[3]==4) # Interval or Left Censored Data
if (length(index)>0){
hfw2$gamma[index]=dwgg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw2$ga[index]=dwga.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw2$gq[index]=dwgq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw2$alpha[index]=dwaa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw2$aq[index]=dwaq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw2$q[index]=dwqq.fun(w2[index],q[index])
hSw2$gamma[index]=Swgg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw2$ga[index]=Swga.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw2$gq[index]=Swgq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw2$alpha[index]=Swaa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw2$aq[index]=Swaq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw2$q[index]=Swqq.fun(w2[index],q[index])
}
}
}
########################################################################
##### Log-likelihood Contribution with Each Obervation with Each Subject
##### and It's First And Second Derivative
########################################################################
if(outtype=="LLF"){
##### Log-likelihood Contribution with Each Obervation with Each Subject
LLC=rep(0,nobs)
## Exact Data
index=which(y[3]==1)
if (length(index)>0) LLC[index]=log(p1[index])+log(fw1[index])-log(alpha[index])-y[index,1]
## Interval Censored Data
index=which(y[3]==3)
if (length(index)>0) LLC[index]=log(p1[index])+log(Sw1[index]-Sw2[index])
## Left Censored Data
index=which(y[3]==4)
if (length(index)>0) LLC[index]=log(p1[index]*(1-Sw2[index]))
## Right Censored Data
index=which(y[3]==2)
if (length(index)>0) LLC[index]=log(p1[index]*Sw1[index]+p2[index])
##### Weighting and Negative Taking of Log-Likelihood Contribution with Each Subject
LLC=-LLC*weight
##### Summation of Log-Likelihood Contribution from All Subjects
TLLC=sum(LLC)
THGNlogLC=TLLC
}else if(outtype=="GLLF"){
##### First Derivative Log-likelihood Contribution with Each Subject
GLLC=list(beta=rep(0,nobs),gamma=rep(0,nobs),alpha=rep(0,nobs),q=rep(0,nobs))
## Exact Data
index=which(y[3]==1)
if (length(index)>0){
base1=fw1[index]
base1g=gfw1$gamma[index]
base1a=gfw1$alpha[index]
base1q=gfw1$q[index]
if (is.reg$beta){ # beta
GLLC$beta[index]=gp1$beta[index]/p1[index]
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1-Aa.fun(x$alpha,par$alpha)[index]/A.fun(x$alpha,par$alpha)[index]
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
## Interval Censored Data
index=which(y[3]==3)
if (length(index)>0){
base1=Sw1[index]-Sw2[index]
base1g=gSw1$gamma[index]-gSw2$gamma[index]
base1a=gSw1$alpha[index]-gSw2$alpha[index]
base1q=gSw1$q[index]-gSw2$q[index]
if (is.reg$beta){ # beta
GLLC$beta[index]=ip1[index]*gp1$beta[index]/p1[index]
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
## Left Censored Data
index=which(y[3]==4)
if (length(index)>0){
base1=p1[index]*(1-Sw2[index])
base1b=ip1[index]*gp1$beta[index]*(1-Sw2[index])
base1g=p1[index]*(-gSw2$gamma[index])
base1a=p1[index]*(-gSw2$alpha[index])
base1q=p1[index]*(-gSw2$q[index])
if (is.reg$beta){ # beta
GLLC$beta[index]=base1b/base1
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
## Right Censored Data
index=which(y[3]==2)
if (length(index)>0){
base1=p1[index]*Sw1[index]+p2[index]
base1b=ip1[index]*gp1$beta[index]*Sw1[index]+ip2[index]*gp2$beta[index]
base1g=p1[index]*gSw1$gamma[index]
base1a=p1[index]*gSw1$alpha[index]
base1q=p1[index]*gSw1$q[index]
if (is.reg$beta){ # beta
GLLC$beta[index]=base1b/base1
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
##### Weighting and Negative Taking of First Derivative Log-Likelihood Contribution with Each Subject
if (is.reg$beta) GLLC$beta=-GLLC$beta*weight*x$beta
GLLC$gamma=-GLLC$gamma*weight*x$gamma
GLLC$alpha=-GLLC$alpha*weight*x$alpha
if(!is.fix.par$q) GLLC$q=-GLLC$q*weight*x$q
##### Summation of First Derivative Log-Likelihood Contribution from All Subjects
TGLLC=NULL
if (is.reg$beta) TGLLC=c(TGLLC,as.vector(apply(GLLC$beta,2,sum))) # beta
TGLLC=c(TGLLC,as.vector(apply(GLLC$gamma,2,sum))) # gamma
TGLLC=c(TGLLC,as.vector(apply(GLLC$alpha,2,sum))) # alpha
if(!is.fix.par$q) TGLLC=c(TGLLC,as.vector(apply(GLLC$q,2,sum))) # q
THGNlogLC=TGLLC
}else if(outtype=="HLLF"){
##### Second Derivative Log-likelihood Contribution with Each Subject
HLLC=list(beta=rep(0,nobs),gamma=rep(0,nobs),alpha=rep(0,nobs),q=rep(0,nobs),
bg=rep(0,nobs),ba=rep(0,nobs),bq=rep(0,nobs),
ga=rep(0,nobs),gq=rep(0,nobs),aq=rep(0,nobs))
## Exact Data
index=which(y[3]==1)
if (length(index)>0){
base1=fw1[index]
base1g=gfw1$gamma[index]
base1a=gfw1$alpha[index]
base1q=gfw1$q[index]
base1gg=hfw1$gamma[index]
base1ga=hfw1$ga[index]
base1gq=hfw1$gq[index]
base1aa=hfw1$alpha[index]
base1aq=hfw1$aq[index]
base1qq=hfw1$q[index]
if (is.reg$beta) HLLC$beta[index]=(hp1$beta[index]*p1[index]-gp1$beta[index]^2)/p1[index]^2
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2-
(Aaa.fun(x$alpha,par$alpha)[index]*A.fun(x$alpha,par$alpha)[index]-Aa.fun(x$alpha,par$alpha)[index]^2)/A.fun(x$alpha,par$alpha)[index]^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
## Right Censored Data
index=which(y[3]==2)
if (length(index)>0){
base1=p1[index]*Sw1[index]+p2[index]
base1b=ip1[index]*gp1$beta[index]*Sw1[index]+ip2[index]*gp2$beta[index]
base1g=p1[index]*gSw1$gamma[index]
base1a=p1[index]*gSw1$alpha[index]
base1q=p1[index]*gSw1$q[index]
base1bb=ip1[index]*hp1$beta[index]*Sw1[index]+ip2[index]*hp2$beta[index]
base1bg=ip1[index]*gp1$beta[index]*gSw1$gamma[index]
base1ba=ip1[index]*gp1$beta[index]*gSw1$alpha[index]
base1bq=ip1[index]*gp1$beta[index]*gSw1$q[index]
base1gg=p1[index]*hSw1$gamma[index]
base1ga=p1[index]*hSw1$ga[index]
base1gq=p1[index]*hSw1$gq[index]
base1aa=p1[index]*hSw1$alpha[index]
base1aq=p1[index]*hSw1$aq[index]
base1qq=p1[index]*hSw1$q[index]
if (is.reg$beta){
HLLC$beta[index]=(base1bb*base1-base1b^2)/base1^2
HLLC$bg[index]=(base1bg*base1-base1b*base1g)/base1^2
HLLC$ba[index]=(base1ba*base1-base1b*base1a)/base1^2
if(!is.fix.par$q) HLLC$bq[index]=(base1bq*base1-base1b*base1q)/base1^2
}
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
## Interval Censored Data
index=which(y[3]==3)
if (length(index)>0){
base1=Sw1[index]-Sw2[index]
base1g=gSw1$gamma[index]-gSw2$gamma[index]
base1a=gSw1$alpha[index]-gSw2$alpha[index]
base1q=gSw1$q[index]-gSw2$q[index]
base1gg=hSw1$gamma[index]-hSw2$gamma[index]
base1ga=hSw1$ga[index]-hSw2$ga[index]
base1gq=hSw1$gq[index]-hSw2$gq[index]
base1aa=hSw1$alpha[index]-hSw2$alpha[index]
base1aq=hSw1$aq[index]-hSw2$aq[index]
base1qq=hSw1$q[index]-hSw2$q[index]
if (is.reg$beta){
HLLC$beta[index]=(hp1$beta[index]*p1[index]-gp1$beta[index]^2)/p1[index]^2
}
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
## Left Censored Data
index=which(y[3]==4)
if (length(index)>0){
base1=p1[index]*(1-Sw2[index])
base1b=ip1[index]*gp1$beta[index]*(1-Sw2[index])
base1g=p1[index]*(-gSw2$gamma[index])
base1a=p1[index]*(-gSw2$alpha[index])
base1q=p1[index]*(-gSw2$q[index])
base1bb=ip1[index]*hp1$beta[index]*(1-Sw2[index])
base1bg=ip1[index]*gp1$beta[index]*(-gSw2$gamma[index])
base1ba=ip1[index]*gp1$beta[index]*(-gSw2$alpha[index])
base1bq=ip1[index]*gp1$beta[index]*(-gSw2$q[index])
base1gg=p1[index]*(-hSw2$gamma[index])
base1ga=p1[index]*(-hSw2$ga[index])
base1gq=p1[index]*(-hSw2$gq[index])
base1aa=p1[index]*(-hSw2$alpha[index])
base1aq=p1[index]*(-hSw2$aq[index])
base1qq=p1[index]*(-hSw2$q[index])
if (is.reg$beta){
HLLC$beta[index]=(base1bb*base1-base1b^2)/base1^2
HLLC$bg[index]=(base1bg*base1-base1b*base1g)/base1^2
HLLC$ba[index]=(base1ba*base1-base1b*base1a)/base1^2
if(!is.fix.par$q) HLLC$bq[index]=(base1bq*base1-base1b*base1q)/base1^2
}
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
##### Weighting and Negative Taking of Second Derivative Log-Likelihood Contribution with Each Subject
if (is.reg$beta){
HLLC$beta=-HLLC$beta*weight*CMultiply(x$beta)
HLLC$bg=-HLLC$bg*weight*CMultiply(x$beta,x$gamma)
HLLC$ba=-HLLC$ba*weight*CMultiply(x$beta,x$alpha)
if(!is.fix.par$q) HLLC$bq=-HLLC$bq*weight*CMultiply(x$beta,x$q)
}
HLLC$gamma=-HLLC$gamma*weight*CMultiply(x$gamma)
HLLC$ga=-HLLC$ga*weight*CMultiply(x$gamma,x$alpha)
if(!is.fix.par$q) HLLC$gq=-HLLC$gq*weight*CMultiply(x$gamma,x$q)
HLLC$alpha=-HLLC$alpha*weight*CMultiply(x$alpha)
if(!is.fix.par$q) HLLC$aq=-HLLC$aq*weight*CMultiply(x$alpha,x$q)
if(!is.fix.par$q) HLLC$q=-HLLC$q*weight*CMultiply(x$q)
##### Summation of Second Derivative Log-Likelihood Contribution from All Subjects
loc2=c(length(par$beta),length(par$beta)+length(par$gamma),length(par$beta)+length(par$gamma)+length(par$alpha),
length(par$beta)+length(par$gamma)+length(par$alpha)+length(par$q))
loc1=loc2-c(length(par$beta),length(par$gamma),length(par$alpha),length(par$q))+1
dim=0
if (is.reg$beta) dim=dim+length(par$beta)
if (is.reg$gamma) dim=dim+length(par$gamma)
if (is.reg$alpha) dim=dim+length(par$alpha)
if (is.reg$q){
if(!is.fix.par$q) dim=dim+length(par$q)
}
THLLC=matrix(0,dim,dim)
### beta
if (is.reg$beta){
# beta
buffer=as.vector(apply(HLLC$beta,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in i:loc2[1]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# bg
buffer=as.vector(apply(HLLC$bg,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in loc1[2]:loc2[2]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# ba
buffer=as.vector(apply(HLLC$ba,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in loc1[3]:loc2[3]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# bq
if(!is.fix.par$q){
buffer=as.vector(apply(HLLC$bq,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in loc1[4]:loc2[4]){
THLLC[i,j]=buffer[k];k=k+1
}
}
}
}
### gamma
buffer=as.vector(apply(HLLC$gamma,2,sum));k=1
for (i in loc1[2]:loc2[2]){ # gamma
for(j in i:loc2[2]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# ga
buffer=as.vector(apply(HLLC$ga,2,sum));k=1
for (i in loc1[2]:loc2[2]){
for(j in loc1[3]:loc2[3]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# gq
if(!is.fix.par$q){
buffer=as.vector(apply(HLLC$gq,2,sum));k=1
for (i in loc1[2]:loc2[2]){
for(j in loc1[4]:loc2[4]){
THLLC[i,j]=buffer[k];k=k+1
}
}
}
### alpha
buffer=as.vector(apply(HLLC$alpha,2,sum));k=1
for (i in loc1[3]:loc2[3]){ # alpha
for(j in i:loc2[3]){
THLLC[i,j]=buffer[k];k=k+1
}
}
if(!is.fix.par$q){ # aq
buffer=as.vector(apply(HLLC$aq,2,sum));k=1
for (i in loc1[3]:loc2[3]){
for(j in loc1[4]:loc2[4]){
THLLC[i,j]=buffer[k];k=k+1
}
}
}
### q
if(!is.fix.par$q){
buffer=as.vector(apply(HLLC$q,2,sum));k=1
for (i in loc1[4]:loc2[4]){ # q
for(j in i:loc2[4]){THLLC[i,j]=buffer[k];k=k+1}}
}
THGNlogLC=t(THLLC)+THLLC;diag(THGNlogLC)=diag(THGNlogLC)/2
}
return(THGNlogLC)
}
#########################
### function: CMultiply()
#########################
CMultiply=function(A,B=NULL){
A=as.matrix(A)
if (is.null(B)){B=A;is.same=T
}else{B=as.matrix(B);is.same=F}
dimA=dim(A);dimB=dim(B)
result=NULL
if (is.same){
for (i in 1:dimA[2]) result=cbind(result,A[,i]*B[,i:dimB[2]])
}else{
if (dimA[1]==dimB[1]){
for (i in 1:dimA[2]) result=cbind(result,A[,i]*B)
}
}
return(result)
}
# CMultiply(x$alpha)
# CMultiply(x$alpha,x$beta)
Try the MixtureRegLTIC package in your browser
Any scripts or data that you put into this service are public.
MixtureRegLTIC documentation built on May 1, 2019, 9:21 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions LLF GLLF HLLF HGNlogL HGNlogLC CMultiply
###################
### function: LLF()
###################
LLF=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype){
LLF=HGNlogL(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="LLF")
return(LLF)
}
####################
### function: GLLF()
####################
GLLF=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype){
GLLF=HGNlogL(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="GLLF")
return(GLLF)
}
####################
### function: HLLF()
####################
HLLF=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype){
HLLF=HGNlogL(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="HLLF")
return(HLLF)
}
#######################
### function: HGNlogL()
#######################
HGNlogL=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="LLF"){
HGNlogL=0
if (length(survtime)==6){
# Calculate truncation contribution
HGNlogL=HGNlogL-HGNlogLC(par.nfix,par.all,is.fix.par,survtime[4:6],x,is.reg,weight,mtype,outtype)
}
# Calculate interval censorship contribution
HGNlogL=HGNlogL+HGNlogLC(par.nfix,par.all,is.fix.par,survtime[1:3],x,is.reg,weight,mtype,outtype)
return(HGNlogL)
}
########################
### function: HGNlogLC()
########################
HGNlogLC=function(par.nfix,par.all,is.fix.par,survtime,x,is.reg,weight,mtype,outtype="LLF"){
##### "par.all" => "par"
par=par.all
##### "par.nfix" => "par"
start=1
if (is.reg$beta){
len=length(par$beta); par$beta=par.nfix[seq(start,start+len-1)]; start=start+len
}
len=length(par$gamma); par$gamma=par.nfix[seq(start,start+len-1)]; start=start+len
len=length(par$alpha); par$alpha=par.nfix[seq(start,start+len-1)]; start=start+len
if (!is.fix.par$q){
len=length(par$q); par$q=par.nfix[seq(start,start+len-1)]; start=start+len
}
##### "y" and "nobs" (Survival Time and Number of Subjects)
y=survtime;nobs=dim(y)[1]
##### Mixture Type (The Type of Distribution Form for Each Subject)
ip2=rep(1,nobs) # "two components with cure"
index=which(mtype==1) # "one component"
if (length(index)>0){ip2[index]=0}
ip1=ip2 # used for derivative with beta
#######################################
##### "alpha" and "q" (scale and shape)
#######################################
alpha=A.fun(x$alpha,par$alpha)
q=Q.fun(x$q,par$q)
##########################################
##### Functions of "p1" & "p2=1-p1"
##### and It's First and Second Derivative
##########################################
##### "p1","p2" (function)
p1=rep(1,nobs)
if(is.reg$beta){
p1=p1.fun(x,par,is.reg)
index=which(mtype==1|mtype==3);if(length(index)>0) p1[index]=1
}
p2=rep(0,nobs);p2=1-p1
if(outtype=="GLLF"|outtype=="HLLF"){
##### "gp1","gp2" (first derivative)
gp1=list(beta=rep(0,nobs))
if (is.reg$beta){
gp1$beta=p1b.fun(x,par,is.reg)
index=which(mtype==1|mtype==3);if(length(index)>0) gp1$beta[index]=0
}
gp2=list(beta=rep(0,nobs));gp2$beta=-gp1$beta
if(outtype=="HLLF"){
##### "hp1", "hp2" (second derivative)
hp1=list(beta=rep(0,nobs))
if (is.reg$beta){
hp1$beta=p1bb.fun(x,par,is.reg)
index=which(mtype==1|mtype==3);if (length(index)>0) hp1$beta[index]=0
}
hp2=list(beta=rep(0,nobs));hp2$beta=-hp1$beta
}
}
################
##### Value of w
################
##### "w1","w2"
w1=rep(NA,nobs);w2=rep(NA,nobs)
index=which(y[,3]==1|y[,3]==3|y[,3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){w1[index]=w.fun(y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)}
index=which(y[,3]==3|y[,3]==4) # Interval or Left Censored Data
if (length(index)>0){w2[index]=w.fun(y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)}
##########################################
##### Density and Survival Functions of w
##### and It's First and Second Derivative
##########################################
##### Density and Survival Functions of w
### "fw1","Sw1","fw2","Sw2"
fw1=rep(NA,nobs);Sw1=rep(NA,nobs)
fw2=rep(NA,nobs);Sw2=rep(NA,nobs)
index=which(y[3]==1|y[3]==3|y[3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){
fw1[index]=dw.fun(w1[index],q[index])
Sw1[index]=Sw.fun(w1[index],q[index])
}
index=which(y[3]==3|y[3]==4) ### Interval Censored or Left Censored Data
if (length(index)>0){
fw2[index]=dw.fun(w2[index],q[index])
Sw2[index]=Sw.fun(w2[index],q[index])
}
if(outtype=="GLLF"|outtype=="HLLF"){
##### First Derivative of Density and Survival Functions of w
### "gfw1","gSw1","gfw2","gSw2"
gfw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
gSw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
gfw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
gSw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs))
index=which(y[3]==1|y[3]==3|y[3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){
gfw1$gamma[index]=dwg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gfw1$alpha[index]=dwa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){
gfw1$q[index]=dwq.fun(w1[index],q[index])
}
gSw1$gamma[index]=Swg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gSw1$alpha[index]=Swa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){gSw1$q[index]=Swq.fun(w1[index],q[index])}
}
index=which(y[3]==3|y[3]==4) # Interval Censored or Left Censored Data
if (length(index)>0){
gfw2$gamma[index]=dwg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gfw2$alpha[index]=dwa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){gfw2$q[index]=dwq.fun(w2[index],q[index])}
gSw2$gamma[index]=Swg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
gSw2$alpha[index]=Swa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q){gSw2$q[index]=Swq.fun(w2[index],q[index])}
}
if(outtype=="HLLF"){
##### Second Derivative of Density and Survival Functions of w
### "hfw1","hSw1","hfw2","hSw2"
hfw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
hSw1=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
hfw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
hSw2=list(gamma=rep(NA,nobs),alpha=rep(NA,nobs),q=rep(NA,nobs),
ga=rep(NA,nobs),gq=rep(NA,nobs),aq=rep(NA,nobs))
index=which(y[3]==1|y[3]==3|y[3]==2) # Exact, Interval Censored or Right Censored Data
if (length(index)>0){
hfw1$gamma[index]=dwgg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw1$ga[index]=dwga.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw1$gq[index]=dwgq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw1$alpha[index]=dwaa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw1$aq[index]=dwaq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw1$q[index]=dwqq.fun(w1[index],q[index])
hSw1$gamma[index]=Swgg.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw1$ga[index]=Swga.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw1$gq[index]=Swgq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw1$alpha[index]=Swaa.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw1$aq[index]=Swaq.fun(w1[index],q[index],y[index,1],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw1$q[index]=Swqq.fun(w1[index],q[index])
}
index=which(y[3]==3|y[3]==4) # Interval or Left Censored Data
if (length(index)>0){
hfw2$gamma[index]=dwgg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw2$ga[index]=dwga.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw2$gq[index]=dwgq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hfw2$alpha[index]=dwaa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw2$aq[index]=dwaq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hfw2$q[index]=dwqq.fun(w2[index],q[index])
hSw2$gamma[index]=Swgg.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw2$ga[index]=Swga.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw2$gq[index]=Swgq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
hSw2$alpha[index]=Swaa.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw2$aq[index]=Swaq.fun(w2[index],q[index],y[index,2],x$gamma[index,],par$gamma,x$alpha[index,],par$alpha)
if(!is.fix.par$q) hSw2$q[index]=Swqq.fun(w2[index],q[index])
}
}
}
########################################################################
##### Log-likelihood Contribution with Each Obervation with Each Subject
##### and It's First And Second Derivative
########################################################################
if(outtype=="LLF"){
##### Log-likelihood Contribution with Each Obervation with Each Subject
LLC=rep(0,nobs)
## Exact Data
index=which(y[3]==1)
if (length(index)>0) LLC[index]=log(p1[index])+log(fw1[index])-log(alpha[index])-y[index,1]
## Interval Censored Data
index=which(y[3]==3)
if (length(index)>0) LLC[index]=log(p1[index])+log(Sw1[index]-Sw2[index])
## Left Censored Data
index=which(y[3]==4)
if (length(index)>0) LLC[index]=log(p1[index]*(1-Sw2[index]))
## Right Censored Data
index=which(y[3]==2)
if (length(index)>0) LLC[index]=log(p1[index]*Sw1[index]+p2[index])
##### Weighting and Negative Taking of Log-Likelihood Contribution with Each Subject
LLC=-LLC*weight
##### Summation of Log-Likelihood Contribution from All Subjects
TLLC=sum(LLC)
THGNlogLC=TLLC
}else if(outtype=="GLLF"){
##### First Derivative Log-likelihood Contribution with Each Subject
GLLC=list(beta=rep(0,nobs),gamma=rep(0,nobs),alpha=rep(0,nobs),q=rep(0,nobs))
## Exact Data
index=which(y[3]==1)
if (length(index)>0){
base1=fw1[index]
base1g=gfw1$gamma[index]
base1a=gfw1$alpha[index]
base1q=gfw1$q[index]
if (is.reg$beta){ # beta
GLLC$beta[index]=gp1$beta[index]/p1[index]
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1-Aa.fun(x$alpha,par$alpha)[index]/A.fun(x$alpha,par$alpha)[index]
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
## Interval Censored Data
index=which(y[3]==3)
if (length(index)>0){
base1=Sw1[index]-Sw2[index]
base1g=gSw1$gamma[index]-gSw2$gamma[index]
base1a=gSw1$alpha[index]-gSw2$alpha[index]
base1q=gSw1$q[index]-gSw2$q[index]
if (is.reg$beta){ # beta
GLLC$beta[index]=ip1[index]*gp1$beta[index]/p1[index]
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
## Left Censored Data
index=which(y[3]==4)
if (length(index)>0){
base1=p1[index]*(1-Sw2[index])
base1b=ip1[index]*gp1$beta[index]*(1-Sw2[index])
base1g=p1[index]*(-gSw2$gamma[index])
base1a=p1[index]*(-gSw2$alpha[index])
base1q=p1[index]*(-gSw2$q[index])
if (is.reg$beta){ # beta
GLLC$beta[index]=base1b/base1
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
## Right Censored Data
index=which(y[3]==2)
if (length(index)>0){
base1=p1[index]*Sw1[index]+p2[index]
base1b=ip1[index]*gp1$beta[index]*Sw1[index]+ip2[index]*gp2$beta[index]
base1g=p1[index]*gSw1$gamma[index]
base1a=p1[index]*gSw1$alpha[index]
base1q=p1[index]*gSw1$q[index]
if (is.reg$beta){ # beta
GLLC$beta[index]=base1b/base1
}
GLLC$gamma[index]=base1g/base1
GLLC$alpha[index]=base1a/base1
if(!is.fix.par$q) GLLC$q[index]=base1q/base1
}
##### Weighting and Negative Taking of First Derivative Log-Likelihood Contribution with Each Subject
if (is.reg$beta) GLLC$beta=-GLLC$beta*weight*x$beta
GLLC$gamma=-GLLC$gamma*weight*x$gamma
GLLC$alpha=-GLLC$alpha*weight*x$alpha
if(!is.fix.par$q) GLLC$q=-GLLC$q*weight*x$q
##### Summation of First Derivative Log-Likelihood Contribution from All Subjects
TGLLC=NULL
if (is.reg$beta) TGLLC=c(TGLLC,as.vector(apply(GLLC$beta,2,sum))) # beta
TGLLC=c(TGLLC,as.vector(apply(GLLC$gamma,2,sum))) # gamma
TGLLC=c(TGLLC,as.vector(apply(GLLC$alpha,2,sum))) # alpha
if(!is.fix.par$q) TGLLC=c(TGLLC,as.vector(apply(GLLC$q,2,sum))) # q
THGNlogLC=TGLLC
}else if(outtype=="HLLF"){
##### Second Derivative Log-likelihood Contribution with Each Subject
HLLC=list(beta=rep(0,nobs),gamma=rep(0,nobs),alpha=rep(0,nobs),q=rep(0,nobs),
bg=rep(0,nobs),ba=rep(0,nobs),bq=rep(0,nobs),
ga=rep(0,nobs),gq=rep(0,nobs),aq=rep(0,nobs))
## Exact Data
index=which(y[3]==1)
if (length(index)>0){
base1=fw1[index]
base1g=gfw1$gamma[index]
base1a=gfw1$alpha[index]
base1q=gfw1$q[index]
base1gg=hfw1$gamma[index]
base1ga=hfw1$ga[index]
base1gq=hfw1$gq[index]
base1aa=hfw1$alpha[index]
base1aq=hfw1$aq[index]
base1qq=hfw1$q[index]
if (is.reg$beta) HLLC$beta[index]=(hp1$beta[index]*p1[index]-gp1$beta[index]^2)/p1[index]^2
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2-
(Aaa.fun(x$alpha,par$alpha)[index]*A.fun(x$alpha,par$alpha)[index]-Aa.fun(x$alpha,par$alpha)[index]^2)/A.fun(x$alpha,par$alpha)[index]^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
## Right Censored Data
index=which(y[3]==2)
if (length(index)>0){
base1=p1[index]*Sw1[index]+p2[index]
base1b=ip1[index]*gp1$beta[index]*Sw1[index]+ip2[index]*gp2$beta[index]
base1g=p1[index]*gSw1$gamma[index]
base1a=p1[index]*gSw1$alpha[index]
base1q=p1[index]*gSw1$q[index]
base1bb=ip1[index]*hp1$beta[index]*Sw1[index]+ip2[index]*hp2$beta[index]
base1bg=ip1[index]*gp1$beta[index]*gSw1$gamma[index]
base1ba=ip1[index]*gp1$beta[index]*gSw1$alpha[index]
base1bq=ip1[index]*gp1$beta[index]*gSw1$q[index]
base1gg=p1[index]*hSw1$gamma[index]
base1ga=p1[index]*hSw1$ga[index]
base1gq=p1[index]*hSw1$gq[index]
base1aa=p1[index]*hSw1$alpha[index]
base1aq=p1[index]*hSw1$aq[index]
base1qq=p1[index]*hSw1$q[index]
if (is.reg$beta){
HLLC$beta[index]=(base1bb*base1-base1b^2)/base1^2
HLLC$bg[index]=(base1bg*base1-base1b*base1g)/base1^2
HLLC$ba[index]=(base1ba*base1-base1b*base1a)/base1^2
if(!is.fix.par$q) HLLC$bq[index]=(base1bq*base1-base1b*base1q)/base1^2
}
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
## Interval Censored Data
index=which(y[3]==3)
if (length(index)>0){
base1=Sw1[index]-Sw2[index]
base1g=gSw1$gamma[index]-gSw2$gamma[index]
base1a=gSw1$alpha[index]-gSw2$alpha[index]
base1q=gSw1$q[index]-gSw2$q[index]
base1gg=hSw1$gamma[index]-hSw2$gamma[index]
base1ga=hSw1$ga[index]-hSw2$ga[index]
base1gq=hSw1$gq[index]-hSw2$gq[index]
base1aa=hSw1$alpha[index]-hSw2$alpha[index]
base1aq=hSw1$aq[index]-hSw2$aq[index]
base1qq=hSw1$q[index]-hSw2$q[index]
if (is.reg$beta){
HLLC$beta[index]=(hp1$beta[index]*p1[index]-gp1$beta[index]^2)/p1[index]^2
}
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
## Left Censored Data
index=which(y[3]==4)
if (length(index)>0){
base1=p1[index]*(1-Sw2[index])
base1b=ip1[index]*gp1$beta[index]*(1-Sw2[index])
base1g=p1[index]*(-gSw2$gamma[index])
base1a=p1[index]*(-gSw2$alpha[index])
base1q=p1[index]*(-gSw2$q[index])
base1bb=ip1[index]*hp1$beta[index]*(1-Sw2[index])
base1bg=ip1[index]*gp1$beta[index]*(-gSw2$gamma[index])
base1ba=ip1[index]*gp1$beta[index]*(-gSw2$alpha[index])
base1bq=ip1[index]*gp1$beta[index]*(-gSw2$q[index])
base1gg=p1[index]*(-hSw2$gamma[index])
base1ga=p1[index]*(-hSw2$ga[index])
base1gq=p1[index]*(-hSw2$gq[index])
base1aa=p1[index]*(-hSw2$alpha[index])
base1aq=p1[index]*(-hSw2$aq[index])
base1qq=p1[index]*(-hSw2$q[index])
if (is.reg$beta){
HLLC$beta[index]=(base1bb*base1-base1b^2)/base1^2
HLLC$bg[index]=(base1bg*base1-base1b*base1g)/base1^2
HLLC$ba[index]=(base1ba*base1-base1b*base1a)/base1^2
if(!is.fix.par$q) HLLC$bq[index]=(base1bq*base1-base1b*base1q)/base1^2
}
HLLC$gamma[index]=(base1gg*base1-base1g^2)/base1^2
HLLC$ga[index]=(base1ga*base1-base1g*base1a)/base1^2
if(!is.fix.par$q) HLLC$gq[index]=(base1gq*base1-base1g*base1q)/base1^2
HLLC$alpha[index]=(base1aa*base1-base1a^2)/base1^2
if(!is.fix.par$q) HLLC$aq[index]=(base1aq*base1-base1a*base1q)/base1^2
if(!is.fix.par$q) HLLC$q[index]=(base1qq*base1-base1q^2)/base1^2
}
##### Weighting and Negative Taking of Second Derivative Log-Likelihood Contribution with Each Subject
if (is.reg$beta){
HLLC$beta=-HLLC$beta*weight*CMultiply(x$beta)
HLLC$bg=-HLLC$bg*weight*CMultiply(x$beta,x$gamma)
HLLC$ba=-HLLC$ba*weight*CMultiply(x$beta,x$alpha)
if(!is.fix.par$q) HLLC$bq=-HLLC$bq*weight*CMultiply(x$beta,x$q)
}
HLLC$gamma=-HLLC$gamma*weight*CMultiply(x$gamma)
HLLC$ga=-HLLC$ga*weight*CMultiply(x$gamma,x$alpha)
if(!is.fix.par$q) HLLC$gq=-HLLC$gq*weight*CMultiply(x$gamma,x$q)
HLLC$alpha=-HLLC$alpha*weight*CMultiply(x$alpha)
if(!is.fix.par$q) HLLC$aq=-HLLC$aq*weight*CMultiply(x$alpha,x$q)
if(!is.fix.par$q) HLLC$q=-HLLC$q*weight*CMultiply(x$q)
##### Summation of Second Derivative Log-Likelihood Contribution from All Subjects
loc2=c(length(par$beta),length(par$beta)+length(par$gamma),length(par$beta)+length(par$gamma)+length(par$alpha),
length(par$beta)+length(par$gamma)+length(par$alpha)+length(par$q))
loc1=loc2-c(length(par$beta),length(par$gamma),length(par$alpha),length(par$q))+1
dim=0
if (is.reg$beta) dim=dim+length(par$beta)
if (is.reg$gamma) dim=dim+length(par$gamma)
if (is.reg$alpha) dim=dim+length(par$alpha)
if (is.reg$q){
if(!is.fix.par$q) dim=dim+length(par$q)
}
THLLC=matrix(0,dim,dim)
### beta
if (is.reg$beta){
# beta
buffer=as.vector(apply(HLLC$beta,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in i:loc2[1]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# bg
buffer=as.vector(apply(HLLC$bg,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in loc1[2]:loc2[2]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# ba
buffer=as.vector(apply(HLLC$ba,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in loc1[3]:loc2[3]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# bq
if(!is.fix.par$q){
buffer=as.vector(apply(HLLC$bq,2,sum));k=1
for (i in loc1[1]:loc2[1]){
for(j in loc1[4]:loc2[4]){
THLLC[i,j]=buffer[k];k=k+1
}
}
}
}
### gamma
buffer=as.vector(apply(HLLC$gamma,2,sum));k=1
for (i in loc1[2]:loc2[2]){ # gamma
for(j in i:loc2[2]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# ga
buffer=as.vector(apply(HLLC$ga,2,sum));k=1
for (i in loc1[2]:loc2[2]){
for(j in loc1[3]:loc2[3]){
THLLC[i,j]=buffer[k];k=k+1
}
}
# gq
if(!is.fix.par$q){
buffer=as.vector(apply(HLLC$gq,2,sum));k=1
for (i in loc1[2]:loc2[2]){
for(j in loc1[4]:loc2[4]){
THLLC[i,j]=buffer[k];k=k+1
}
}
}
### alpha
buffer=as.vector(apply(HLLC$alpha,2,sum));k=1
for (i in loc1[3]:loc2[3]){ # alpha
for(j in i:loc2[3]){
THLLC[i,j]=buffer[k];k=k+1
}
}
if(!is.fix.par$q){ # aq
buffer=as.vector(apply(HLLC$aq,2,sum));k=1
for (i in loc1[3]:loc2[3]){
for(j in loc1[4]:loc2[4]){
THLLC[i,j]=buffer[k];k=k+1
}
}
}
### q
if(!is.fix.par$q){
buffer=as.vector(apply(HLLC$q,2,sum));k=1
for (i in loc1[4]:loc2[4]){ # q
for(j in i:loc2[4]){THLLC[i,j]=buffer[k];k=k+1}}
}
THGNlogLC=t(THLLC)+THLLC;diag(THGNlogLC)=diag(THGNlogLC)/2
}
return(THGNlogLC)
}
#########################
### function: CMultiply()
#########################
CMultiply=function(A,B=NULL){
A=as.matrix(A)
if (is.null(B)){B=A;is.same=T
}else{B=as.matrix(B);is.same=F}
dimA=dim(A);dimB=dim(B)
result=NULL
if (is.same){
for (i in 1:dimA[2]) result=cbind(result,A[,i]*B[,i:dimB[2]])
}else{
if (dimA[1]==dimB[1]){
for (i in 1:dimA[2]) result=cbind(result,A[,i]*B)
}
}
return(result)
}
# CMultiply(x$alpha)
# CMultiply(x$alpha,x$beta)
Try the MixtureRegLTIC package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.