R/EM.PScr.R

In PScr: Estimation for the Power Series Cure Rate Model

EM.PScr <-
function(t, delta, z, model=1, dist=1, max.iter=1000, prec=1e-4)
{
passoE<-function(psi,data,model=1,dist=1,q=1)
{
t<-data$t
delta<-data$delta
z<-matrix(data$z,nrow=r)
beta<-matrix(psi[1:r],ncol=1)
alpha<-psi[r+1]
sigma<-psi[r+2]
if(model==3 || model==5)
{
q<-psi[r+3]
}
if(dist==1)
{
F<-2*pnorm(t/sigma)-1-exp((1/2)*(alpha*log(2)-log(pi))+alpha*(log(sigma)-log(t))+lgamma((alpha+1)/2)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE))
}
if(dist==2)
{
F<-1-exp(-(t/sigma)^alpha)
}
if(dist==3)
{
F<-pgamma(t,shape=alpha,scale=sigma)
}
if(dist==4)
{
F<-pbisa(t,shape=alpha,scale=sigma)
}
if(dist==5)
{
F<-plnorm(t,meanlog=alpha,sdlog=sigma)
}
S<-1-F
if(model==1 || model==4)
{
theta=exp(t(z)%*%beta)
}
if(model==2 || model==3 || model==5 || model==6)
{
theta=plogis(t(z)%*%beta)
}
M<-c()
if(model==1)
{
M=delta+theta*S
}
if(model==2)
{
mu=theta*S
num=(1-mu)^2*(log(1-mu))^2-mu*(2*mu+delta-2)*log(1-mu)+(1-delta)*mu^2
den=-(1-mu)*log(1-mu)*((1-mu)*log(1-mu)+mu)
M=num/den
}
if(model==3)
{
mu=theta*S
M=(delta+q*mu)/(1-mu)
}
if(model==4)
{
mu=theta*S
M=(delta+q*mu)/(1+mu)
}
if(model==5)
{
mu1=polyloga(theta*S,q-1)/polyloga(theta*S,q)-1
mu2=polyloga(theta*S,q-2)/polyloga(theta*S,q)-2*polyloga(theta*S,q-1)/polyloga(theta*S,q)+1
M=(1-delta)*mu1+delta*mu2/mu1
}
if(model==6)
{
mu=theta*S
M=(1-delta)*2*mu*(1-mu)^(-1)+delta*(3*mu*(1-mu)^(-1)+1)
}
return(M=c(M))
}
cura.inicial=function(p,model=1,q=1)
{
beta=c()
if(model==1)
{
beta[1]=log(p[1])
if(length(p)>1)
{
for(i in 2:length(p))
{
beta[i]=log(p[i])-sum(beta[1:(i-1)])
}
}
}
if(model==2)
{
cure.log<-function(x,p)
{
theta=exp(x)/(1+exp(x))
-theta/log(1-theta)-p
}
beta=try(uniroot(cure.log, c(-10, 10),p=p[1]))$root
if(length(p)>1)
{
for(i in 2:length(p))
{
beta[i]=try(uniroot(cure.log, c(-10, 10),p=p[i]))$root-sum(beta[1:(i-1)])
}
}
}
if(model==3)
{
beta[1]=log(1-p[1]^(1/q))
if(length(p)>1)
{
for(i in 2:length(p))
{
beta[i]=log(1-p[i]^(1/q))-sum(beta[1:(i-1)])
}
}
}
if(model==4)
{
beta[1]=log(p[1]^(-1/q)-1)
if(length(p)>1)
{
for(i in 2:length(p))
{
beta[i]=log(p[i]^(-1/q)-1)-sum(beta[1:(i-1)])
}
}
}
return(beta)
}
Q1<-function(beta,M,data,model=1,dist=1,q=1)
{
t<-data$t
delta<-data$delta
z<-matrix(data$z,nrow=r)
beta<-matrix(beta,ncol=1)
if(model==1 || model==4)
{
theta=exp(t(z)%*%beta)
}
if(model==2 || model==3 || model==5 || model==6)
{
theta=plogis(t(z)%*%beta)
}
if(model==1)
{
log.A=theta
}
if(model==2)
{
log.A=log(-log(1-theta)/theta)
}
if(model==3)
{
log.A=-q*log(1-theta)
}
if(model==4)
{
log.A=q*log(1+theta)
}
if(model==5)
{
log.A=log(polyloga(theta,q))-log(theta)
}
if(model==6)
{
log.A=-2*log1p(-theta)
}
-sum(M*log(theta)-log.A)
}
Q2<-function(loglambda,M,data,model=1,dist=1)
{
t<-data$t
delta<-data$delta
z<-matrix(data$z,nrow=r)
if(dist==2 || dist==5)
{
alpha=loglambda[1]
sigma=exp(loglambda[2])
}
if(dist==1 || dist==3 || dist==4)
{
lambda=exp(loglambda)
alpha=lambda[1]
sigma=lambda[2]
}
if(dist==1)
{
F<-2*pnorm(t/sigma)-1-exp((1/2)*(alpha*log(2)-log(pi))+alpha*(log(sigma)-log(t))+lgamma((alpha+1)/2)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE))
logf<-log(alpha)+(1/2)*(alpha*log(2)-log(pi))+alpha*log(sigma)+lgamma((alpha+1)/2)-(alpha+1)*log(t)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE)
logS<-log1p(-F)
}
if(dist==2)
{
logS<--(t/sigma)^alpha
logf<-log(alpha)-log(sigma)+(alpha-1)*(log(t)-log(sigma))-(t/sigma)^alpha
}
if(dist==3)
{
logS<-pgamma(t,shape=alpha,scale=sigma,log.p=TRUE,lower.tail=FALSE)
logf<-dgamma(t,shape=alpha,scale=sigma,log=TRUE)
}
if(dist==4)
{
logS<-pbisa(t,shape=alpha,scale=sigma,log.p=TRUE,lower.tail=FALSE)
logf<-dbisa(t,shape=alpha,scale=sigma,log=TRUE)
}
if(dist==5)
{
logS<-plnorm(t,meanlog=alpha,sdlog=sigma,log.p=TRUE,lower.tail=FALSE)
logf<-dlnorm(t,meanlog=alpha,sdlog=sigma,log=TRUE)
}
-sum((M-delta)*logS+delta*logf)
}
Q3<-function(logq,psi,data,model,dist)
{
q=exp(logq)
t<-data$t
delta<-data$delta
z<-matrix(data$z,nrow=r)
beta<-matrix(psi[1:r],ncol=1)
alpha<-psi[r+1]
sigma<-psi[r+2]
if(dist==1)
{
F<-2*pnorm(t/sigma)-1-exp((1/2)*(alpha*log(2)-log(pi))+alpha*(log(sigma)-log(t))+lgamma((alpha+1)/2)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE))
logf<-log(alpha)+(1/2)*(alpha*log(2)-log(pi))+alpha*log(sigma)+lgamma((alpha+1)/2)-(alpha+1)*log(t)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE)
}
if(dist==2)
{
F<-1-exp(-(t/sigma)^alpha)
logf<-log(alpha)-log(sigma)+(alpha-1)*(log(t)-log(sigma))-(t/sigma)^alpha
}
if(dist==3)
{
F<-pgamma(t,shape=alpha,scale=sigma)
logf<-dgamma(t,shape=alpha,scale=sigma,log=TRUE)
}
if(dist==4)
{
F<-pbisa(t,shape=alpha,scale=sigma)
logf<-dbisa(t,shape=alpha,scale=sigma,log=TRUE)
}
if(dist==5)
{
F<-plnorm(t,meanlog=alpha,sdlog=sigma)
logf<-dlnorm(t,meanlog=alpha,sdlog=sigma,log=TRUE)
}
S<-1-F
if(model==4)
{
theta=exp(t(z)%*%beta)
}
if(model==3 || model==5 || model==6)
{
theta=plogis(t(z)%*%beta)
}
if(model==3)
{
mu=theta*S
log.A=-q*log1p(-theta)
log.AS=-q*log1p(-mu)
log.fpop=log(q)+(q-1)*(log1p(-theta)-log1p(-theta*S))
}
if(model==4)
{
mu=theta*S
log.A=q*log(1+theta)
log.AS=q*log(1+mu)
log.fpop=(q-1)*log(1/(1+theta)+theta*S/(1+theta))+log(theta)-log(1+theta)+logf
}
if(model==5)
{
mu=theta*S
log.A=log(polyloga(theta,q))
log.AS=log(polyloga(mu,q))
log.fpop=-log(polyloga(theta,q))+log(polyloga(mu,q-1)-polyloga(mu,q))
}
if(model==6)
{
mu=theta*S
log.A=-2*log1p(-theta)
log.AS=-2*log1p(-mu)
log.fpop=log(2)+2*log1p(-theta)+log(theta)+logf-3*log1p(-mu)
}
log.Spop=log.AS-log.A
-sum(delta*log.fpop+(1-delta)*log.Spop)
}
buscar.indice<-function(vector,a)
{
ind=c()
for(i in 1:length(vector))
{
if(vector[i]==a)
{ind=c(ind,i)}
}
ind
}
llikeobserved<-function(psi,data,model,dist)
{
t<-data$t
delta<-data$delta
z<-matrix(data$z,nrow=r)
beta<-matrix(psi[1:r],ncol=1)
alpha<-psi[r+1]
sigma<-psi[r+2]
if(model==3 || model==5)
{
q=psi[r+3]
}
if(dist==1)
{
F<-2*pnorm(t/sigma)-1-exp((1/2)*(alpha*log(2)-log(pi))+alpha*(log(sigma)-log(t))+lgamma((alpha+1)/2)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE))
logf<-log(alpha)+(1/2)*(alpha*log(2)-log(pi))+alpha*log(sigma)+lgamma((alpha+1)/2)-(alpha+1)*log(t)+pgamma(t^2/(2*sigma^2),shape=(alpha+1)/2,rate=1,log.p=TRUE)
}
if(dist==2)
{
F<-1-exp(-(t/sigma)^alpha)
logf<-log(alpha)-log(sigma)+(alpha-1)*(log(t)-log(sigma))-(t/sigma)^alpha
}
if(dist==3)
{
F<-pgamma(t,shape=alpha,scale=sigma)
logf<-dgamma(t,shape=alpha,scale=sigma,log=TRUE)
}
if(dist==4)
{
F<-pbisa(t,shape=alpha,scale=sigma)
logf<-dbisa(t,shape=alpha,scale=sigma,log=TRUE)
}
if(dist==5)
{
F<-plnorm(t,meanlog=alpha,sdlog=sigma)
logf<-dlnorm(t,meanlog=alpha,sdlog=sigma,log=TRUE)
}
S<-1-F
if(model==1 || model==4)
{
theta=exp(t(z)%*%beta)
}
if(model==2 || model==3 || model==5 || model==6)
{
theta=plogis(t(z)%*%beta)
}
if(model==1)
{
mu=theta*S
log.A=theta
log.AS=mu
log.fpop=log(theta)+logf-theta*F
}
if(model==2)
{
mu=theta*S
log.A=log(-log(1-theta))-log(theta)
log.AS=log(-log(1-mu))-log(mu)
log.fpop=logf+log(-1/(S*log(1-theta)))+log(log(1-theta*S)/S+theta/(1-theta*S))
}
if(model==3)
{
mu=theta*S
log.A=-q*log1p(-theta)
log.AS=-q*log1p(-mu)
log.fpop=log(q)+(q-1)*(log1p(-theta)-log1p(-theta*S))+log(theta)+log1p(-theta)+logf-2*log1p(-theta*S)
}
if(model==4)
{
mu=theta*S
log.A=q*log(1+theta)
log.AS=q*log(1+mu)
log.fpop=(q-1)*log(1/(1+theta)+theta*S/(1+theta))+log(theta)-log(1+theta)+logf
}
if(model==5)
{
mu=theta*S
log.A=log(polyloga(theta,q))+log(S)
log.AS=log(polyloga(mu,q))
logh=logf-log(S)
log.fpop=logh-log(S)-log(polyloga(theta,q))+log(polyloga(mu,q-1)-polyloga(mu,q))
}
if(model==6)
{
mu=theta*S
log.A=-2*log1p(-theta)
log.AS=-2*log1p(-mu)
logh=logf-log(S)
log.fpop=log(2)+2*log1p(-theta)+log(theta)+logf-3*log1p(-mu)
}
log.Spop=log.AS-log.A
-sum(delta*log.fpop+(1-delta)*log.Spop)
}
polyloga<-function(theta,q)
{
theta*lerch(x=theta, s=q, v=1,tolerance = 1e-10,iter=100)
}
q=1
r=min(dim(z))
data=list(t=t,delta=delta,z=z)
KM=survfit(Surv(t,delta)~1)$surv
KM=KM[length(KM)]
psi=c(qlogis(KM),rep(0,nrow(z)-1),1,sqrt(mean(t^2/2)))
if(model==3 || model==5)
{
psi=c(psi, q)
}
i=1; known=0
while(i<=max.iter)
{
M=passoE(psi,data,model,dist)
beta=psi[1:r]
lambda=psi[(r+1):(r+2)]
if(model==3 || model==5)
{
q=psi[r+3]
}
maximo1=nlminb(beta,Q1,M=M,data=data,model=model,dist=dist,q=q)
beta=maximo1$par
if(dist==2 || dist==5)
{
maximo2=nlminb(c(lambda[1],log(lambda[2])),Q2,M=M,data=data,model=model,dist=dist)
lambda=c(maximo2$par[1],exp(maximo2$par)[2])
}
if(dist==1 || dist==3 || dist==4)
{
maximo2=nlminb(log(lambda),Q2,M=M,data=data,model=model,dist=dist)
lambda=exp(maximo2$par)
}
if((model==3  || model==5) && known==0)
{
maximo3=nlminb(log(q),Q3,data=data,model=model,dist=dist,psi=c(beta,lambda))
q=exp(maximo3$par)
}
dif=-llikeobserved(c(beta,lambda,q),data,model,dist)+llikeobserved(psi,data,model,dist)
if(dif>prec)
{
psi=c(beta,lambda)
if(model==3  || model==5)
{
psi=c(beta,lambda,q)
}
}
if(dif<prec)
{
i=max.iter
}
i=i+1
}
if(model!=3)
{
varianza.psi=solve(hessian(llikeobserved,x0=psi,data=data,model=model,dist=dist))
}
if(model==3 || model==5)
{
varianza.psi=solve(hessian(llikeobserved,x0=psi,data=data,model=model,dist=dist))
}
estimates<-matrix(c(psi,sqrt(diag(varianza.psi))),ncol=2)
colnames(estimates)=c("Estimate","S.e.")
if(model==1 || model==2 || model==4)
{
rownames(estimates)=c(paste("beta",0:(r-1)),"alpha","sigma")
}
if(model==3 || model==5)
{
rownames(estimates)=c(paste("beta",0:(r-1)),"alpha","sigma","q")
}
loglike=llikeobserved(psi, data, model, dist)
list(estimate=estimates, loglike=loglike, AIC=2*loglike+2*length(psi), 
BIC=2*loglike+length(psi)*log(length(t)))
}