#This program will fit the conventional Generalized Gamma regression model and return the value of 
#estimated betas, sigma, lambda, and -2 log likelihood

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 #

data.0<-full.data[full.data[,4]==0, ] #all censored times
data.1<-full.data[full.data[,4]==1, ] #all uncensored times

m2loglikgralgammaconventional<-function(para)
{
beta.1 <- para[1]*data.1[,5]+para[2]*data.1[,6]+para[3]*data.1[,7]+para[4]*data.1[,8]
beta.0 <- para[1]*data.0[,5]+para[2]*data.0[,6]+para[3]*data.0[,7]+para[4]*data.0[,8]
sigma <- para[5]
#for general gamma
lambda <- para[6]

numuncens <- nrow(data.1)
numcens <- nrow(data.0)

#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.1,sigma,lambda,data.1[,3])/sggamma(beta.1,sigma,lambda,data.1[,2])))
#censored event
loglik<-loglik+sum(log(sggamma(beta.0,sigma,lambda,data.0[,3])/sggamma(beta.0,sigma,lambda,data.0[,2])))
(-2)*loglik
}

#for general gamma
ggconv.mle_nlminb(start=c(0.5,0.5,2,3,1,1), obj=m2loglikgralgammaconventional)
print(c("GENERAL GAMMA, CONVENTIONAL REGRESSION"))
print(c("beta=",round(ggconv.mle$para[1:4],3)))
print(c("sigma=",round(ggconv.mle$para[5],3)))
#for general gamma
print(c("lambda=",round(ggconv.mle$para[6],3)))

#-2 log likelihood includes the negative of the sum of the logarithms of the uncensored times to event
print(c("-2 Log Likelihood=",round(ggconv.mle$objective,2)))