More than 3 years have passed since last update.

EVENT HISTORY ANALYSIS IN DEMOGRAPHY, DANIEL et.al, CLARENDON Press

RStudio

Posted at 2021-09-05


# Method of Estimation of Parametric Models p.145

# Exponential model

x1=c(rep(0,10),1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0)

x2=c(2,3,0,3,2,1,1,1,3,0,2,0,2,0,2,3,1,2,2,1,1,2,1,3,3,1,1,2,3,1,1,3,3,3,1,3,0,1,0,2,3,2,0,2,1,2,0,3,1,1,3,3,3,2,0,3,3,2,3,3,0,3,0,2,2,2,3,1,3,0,0,3,3,3,0,2,0,1,3,0)

x3=c(1,0,1,1,1,1,0,1,0,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,0,1)

t=c(15,13,15,10,3,2,2,2,1,18,4,5,2,6,1,1,22,11,19,1,23,1,22,22,10,1,31,1,3,4,10,20,2,6,9,1,2,1,3,8,31,1,30,2,5,2,6,5,5,10,3,7,1,1,5,1,2,11,21,14,3,1,2,6,31,4,14,25,6,31,11,7,20,3,8,13,8,1,14,31)

delta=c(rep(1,9),0,rep(1,6),0,0,0,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,1,0,rep(1,9),0,1,0,1,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,0)

times=sort(unique(t))

Z=cbind(x1,x2,x3)


ite=120

eta=1

beta=rep(1,ncol(X))

for(l in 1:ite){
  
dbeta=apply(Z[delta==1,],2,sum)-t(Z)%*%c(t*exp(Z%*%beta))

ddbeta=t(Z)%*%diag(c(t*exp(Z%*%beta)))%*%Z

beta=beta+eta*solve(ddbeta)%*%dbeta

loglik=sum(Z%*%beta-t*exp(Z%*%beta))

print(loglik)
  
}


# Method of Estimation of Parametric Models p.145

# The Weibull model

x1=c(rep(0,10),1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0)

x2=c(2,3,0,3,2,1,1,1,3,0,2,0,2,0,2,3,1,2,2,1,1,2,1,3,3,1,1,2,3,1,1,3,3,3,1,3,0,1,0,2,3,2,0,2,1,2,0,3,1,1,3,3,3,2,0,3,3,2,3,3,0,3,0,2,2,2,3,1,3,0,0,3,3,3,0,2,0,1,3,0)

x3=c(1,0,1,1,1,1,0,1,0,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,0,1)

t=c(15,13,15,10,3,2,2,2,1,18,4,5,2,6,1,1,22,11,19,1,23,1,22,22,10,1,31,1,3,4,10,20,2,6,9,1,2,1,3,8,31,1,30,2,5,2,6,5,5,10,3,7,1,1,5,1,2,11,21,14,3,1,2,6,31,4,14,25,6,31,11,7,20,3,8,13,8,1,14,31)

delta=c(rep(1,9),0,rep(1,6),0,0,0,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,1,0,rep(1,9),0,1,0,1,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,0)

times=sort(unique(t))

Z=cbind(x1,x2,x3)


ite=10000

eta=0.00001

beta=rep(1,ncol(X))

p=1

lam=1

for(l in 1:ite){
  
dbeta=apply(Z[delta==1,],2,sum)-t(Z)%*%c(((p*t)^lam)*exp(Z%*%beta))

dp=lam*sum(delta)/p-lam*(p^(lam-1))*sum(exp(Z%*%beta)*t^lam)

dlam=sum(delta)/lam+sum(delta)*log(p)+sum(log(t[delta==1]))-(p^lam)*sum((t^lam)*log(p*t)*exp(Z%*%beta))

df=c(dlam,dp,dbeta)

param=c(lam,p,beta)+eta*df

lam=param[1];p=param[2];beta=param[-c(1:2)]

loglik=sum((log(lam)+log(p)+(lam-1)*log(p*t)+Z%*%beta)[delta==1])-sum(exp(Z%*%beta)*(p*t)^lam)

print(loglik)
  
}


# Method of Estimation of Parametric Models p.145

# The Gompertz model

x1=c(rep(0,10),1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0)

x2=c(2,3,0,3,2,1,1,1,3,0,2,0,2,0,2,3,1,2,2,1,1,2,1,3,3,1,1,2,3,1,1,3,3,3,1,3,0,1,0,2,3,2,0,2,1,2,0,3,1,1,3,3,3,2,0,3,3,2,3,3,0,3,0,2,2,2,3,1,3,0,0,3,3,3,0,2,0,1,3,0)

x3=c(1,0,1,1,1,1,0,1,0,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,0,1)

t=c(15,13,15,10,3,2,2,2,1,18,4,5,2,6,1,1,22,11,19,1,23,1,22,22,10,1,31,1,3,4,10,20,2,6,9,1,2,1,3,8,31,1,30,2,5,2,6,5,5,10,3,7,1,1,5,1,2,11,21,14,3,1,2,6,31,4,14,25,6,31,11,7,20,3,8,13,8,1,14,31)

delta=c(rep(1,9),0,rep(1,6),0,0,0,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,1,0,rep(1,9),0,1,0,1,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,0)

times=sort(unique(t))

Z=cbind(x1,x2,x3)


ite=30000

eta=10^(-5)

beta=rep(0.01,ncol(X))

p=0.1

lam=0.1

for(l in 1:ite){
  
dbeta=apply(Z[delta==1,],2,sum)+t(Z)%*%c((1-exp(p*t))*exp(Z%*%beta))*lam

dp=sum(delta)/p+sum(t[delta==1])-lam*sum(t*exp(p*t)*exp(Z%*%beta))

dlam=sum(delta)/lam+sum((1-exp(p*t))*exp(Z%*%beta))

df=c(dlam,dp,dbeta)

param=c(lam,p,beta)+eta*df

lam=param[1];p=param[2];beta=param[-c(1:2)]

loglik=sum(exp(Z%*%beta)*(1-exp(p*t)))*lam+sum((log(lam)+log(p)+(p*t)+Z%*%beta)[delta==1])

print(loglik)
  
}


# Method of Estimation of Parametric Models p.145

# Log-logistic model

x1=c(rep(0,10),1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0)

x2=c(2,3,0,3,2,1,1,1,3,0,2,0,2,0,2,3,1,2,2,1,1,2,1,3,3,1,1,2,3,1,1,3,3,3,1,3,0,1,0,2,3,2,0,2,1,2,0,3,1,1,3,3,3,2,0,3,3,2,3,3,0,3,0,2,2,2,3,1,3,0,0,3,3,3,0,2,0,1,3,0)

x3=c(1,0,1,1,1,1,0,1,0,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,0,1)

t=c(15,13,15,10,3,2,2,2,1,18,4,5,2,6,1,1,22,11,19,1,23,1,22,22,10,1,31,1,3,4,10,20,2,6,9,1,2,1,3,8,31,1,30,2,5,2,6,5,5,10,3,7,1,1,5,1,2,11,21,14,3,1,2,6,31,4,14,25,6,31,11,7,20,3,8,13,8,1,14,31)

delta=c(rep(1,9),0,rep(1,6),0,0,0,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,1,0,rep(1,9),0,1,0,1,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,0)

times=sort(unique(t))

Z=cbind(x1,x2,x3)


ite=3000

eta=10^(-3)

beta=rep(0.01,ncol(X))

p=0.1

lam=0.1

for(l in 1:ite){
  
dbeta=t(Z[delta==1,])%*%c(1/(1+exp(Z%*%beta)*t^lam)[delta==1])-t(Z)%*%c(1/(1+exp(-Z%*%beta)*t^(-lam)))

dlam=sum(delta)/lam+sum((log(t)/(1+exp(Z%*%beta)*t^lam))[delta==1])-sum((log(t)/(1+exp(-Z%*%beta)*t^(-lam))))

df=c(dlam,dbeta)

param=c(lam,beta)+eta*df

lam=param[1];beta=param[-c(1)]

loglik=-sum(log(1+exp(Z%*%beta)*t^(lam)))+sum((log(lam)-log(t)-log(1+exp(-Z%*%beta)*t^(-lam)))[delta==1])

print(loglik)
  
}

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up