#This program will fit the saturated Generalized Gamma regression model and return the value of 
#estimated betas, estimated sigmas, estimated lambdas, and -2 log likelihood for each period
#Also, the -2 log likelihood for each period will be summed over all four periods to get a final value for the
#saturated model

full.data<-matrix(scan("E:/GeneralGamma/fulldata.dat"), ncol=9,byrow=T)
#fulldata.dat is file with full data on nine variables for 2314 records, need to update path
#9 variables: publicID entry exit event p1 p2 p3 p4 period1234 #

m2loglikgralgamma<-function(para)
{
beta <- para[1]
sigma <- para[2]
lambda <- para[3]

#Density Function for the Generalized Gamma Distribution
dggamma <- function(beta, sigma, lambda, timet)
{
s <- lambda^(-2)
r <- s * exp(( - beta * lambda)/sigma)

# consider timet below exp(-500) at exp(-500) to avoid division by zero
timet[timet<exp(-500)] <- exp(-500)


b <- lambda/sigma
dens <- (lambda>0)*((b*(timet^(b-1))*dgamma(timet^b, shape=s, rate=r)))
dens <- dens+(lambda<0)*((-b*(timet^(b-1))*dgamma(timet^b, shape=s, rate=r)))
dens
}

#Survival Function for the Generalized Gamma Distribution
sggamma <- function(beta, sigma, lambda, timet)
{
s <- lambda^(-2)
r <- s * exp(( - beta * lambda)/sigma)

# consider timet below exp(-500) at exp(-500) to avoid division by zero
timet[timet<exp(-500)] <- exp(-500)


b <- lambda/sigma
surv <- (lambda>0)*(1-pgamma(timet^b, shape=s, rate=r))
surv <- surv+(lambda<0)*pgamma(timet^b, shape=s, rate=r)
surv[surv<exp(-500)] <- exp(-500)
surv
}

loglik<-0
#uncensored
loglik<-loglik+sum(log(dggamma(beta,sigma,lambda,data.1[,3])/sggamma(beta,sigma,lambda,data.1[,2])))
#censored event
loglik<-loglik+sum(log(sggamma(beta,sigma,lambda,data.0[,3])/sggamma(beta,sigma,lambda,data.0[,2])))
(-2)*loglik
}

#Period 1
data<-full.data[full.data[,5]==1,] #select period 1, entry and exit time, event;
data.0<-data[data[,4]==0,] #censored times
data.1<-data[data[,4]==1,] #uncensored times

#minimize (-2) log likelihood starting at start=
gg.mle_nlminb(start=c(0,1,0.5), obj=m2loglikgralgamma)
print(c("GENERAL GAMMA FOR PERIOD 1"))
print(c("beta=",round(gg.mle$para[1],3)))
print(c("sigma=",round(gg.mle$para[2],3)))
print(c("lambda=",round(gg.mle$para[3],3)))
print(c("-2 Log Likelihood=",round(gg.mle$objective,2)))
totalm2loglik_gg.mle$objective

#Period 2
data<-full.data[full.data[,6]==1,] #select period 2, entry and exit time, event;
data.0<-data[data[,4]==0,] #censored times
data.1<-data[data[,4]==1,] #uncensored times

#minimize (-2) log likelihood starting at start=
gg.mle_nlminb(start=c(0,1,0.5), obj=m2loglikgralgamma)
print(c("GENERAL GAMMA FOR PERIOD 2"))
print(c("beta=",round(gg.mle$para[1],3)))
print(c("sigma=",round(gg.mle$para[2],3)))
print(c("lambda=",round(gg.mle$para[3],3)))
print(c("-2 Log Likelihood=",round(gg.mle$objective,2)))
totalm2loglik_totalm2loglik + gg.mle$objective

#Period 3
data<-full.data[full.data[,7]==1,] #select period 3, entry and exit time, event;
data.0<-data[data[,4]==0, ] #censored times
data.1<-data[data[,4]==1, ] #uncensored times

#minimize (-2) log likelihood starting at start=
gg.mle_nlminb(start=c(0,1,0.5), obj=m2loglikgralgamma)
print(c("GENERAL GAMMA FOR PERIOD 3"))
print(c("beta=",round(gg.mle$para[1],3)))
print(c("sigma=",round(gg.mle$para[2],3)))
print(c("lambda=",round(gg.mle$para[3],3)))
print(c("-2 Log Likelihood=",round(gg.mle$objective,2)))
totalm2loglik_totalm2loglik + gg.mle$objective

#Period 4
data<-full.data[full.data[,8]==1,] #select period 4, entry and exit time, event;
data.0<-data[data[,4]==0, ] #censored times
data.1<-data[data[,4]==1, ] #uncensored times

#minimize (-2) log likelihood starting at start=
gg.mle_nlminb(start=c(0,1,0.5), obj=m2loglikgralgamma)
print(c("GENERAL GAMMA FOR PERIOD 4"))
print(c("beta=",round(gg.mle$para[1],3)))
print(c("sigma=",round(gg.mle$para[2],3)))
print(c("lambda=",round(gg.mle$para[3],3)))
print(c("-2 Log Likelihood=",round(gg.mle$objective,2)))
totalm2loglik_totalm2loglik + gg.mle$objective
#-2 log likelihood includes the negative of the sum of the logarithms of the uncensored times to event
print(c("Total -2 Log Likelihood=",round(totalm2loglik,2)))