#Program to analyze the data in Table 11.1 for Example 11.1
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table111.txt")
count=yd[,1]
x1 = yd[,2]
x2 = yd[,3]
x3 = yd[,4]
x4 = yd[,5]
n = length(x1)
x0= rep(1, n)
#Poisson regression model
md1 = count~x1 + x2 + x3 + x4
gl1 = glm(md1, family=poisson("log"))
anova(gl1, test="Chisq")
summary(gl1)

xx = cbind(x0,x1,x2,x3,x4)         
#To find the inital estimates for the regression coefficients
par = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(count+0.5))

#Fitting Poisson regression model
#f.po is the function that defines the (negative) log likelihood 
#This log likelihood includes the constant part with no parameters
f.po <- function(para.p, xx, count)
{  
m1.x = xx%*%para.p
mu.x = exp(m1.x)
-sum(count*log(mu.x) - mu.x - lgamma(count+1))
}
#The function nlminb does not return the hessian matrix
I.p = nlminb(start=par,objective=f.po,xx=xx,count=count)
print("Poisson Regression Model")
I.p[[1]] #parameter estimates
I.p[[2]] #minus log-likelihood function
I.p[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.op=nlm(f.po,par,xx=xx,count=count,hessian=TRUE)
est=I.op$estimate
log.l=I.op$minimum
hess=I.op$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "x1","x2","x3","x4")
log.l
res

#Poisson regression model
md2 = count~x1 + x3 + x4
gl2 = glm(md2, family=poisson("log"))
anova(gl2, test="Chisq")
summary(gl2)

##########################################################################

#Program to analyze the data in Table 11.2 for Example 11.2
yd=read.table("C:\\Research\\Book\\Data\\table112.txt")
age = yd[,1]
alc = yd[,2]
yes = yd[,3]
noo = yd[,4]
tot = yes + noo
age2=age^2
age.f=as.factor(age)
alc.f=as.factor(alc)
n.alc = 2-as.numeric(alc)

#Poisson regression model
md1 = yes~age.f + alc.f
gl1 = glm(md1, family=poisson("log"), offset=log(tot))
anova(gl1, test="Chisq")
summary(gl1)

#Poisson regression model
md2 = yes~age + age2 + n.alc
gl2 = glm(md2, family=poisson("log"), offset=log(tot))
anova(gl2, test="Chisq")
summary(gl2)
pred=predict(gl2,type="response")
res.dev = residuals(gl2,type="deviance")
res.chi = residuals(gl2,type="pearson")

cbind(age,alc,yes,pred,res.chi,res.dev)

#Logistic regression model
yy=cbind(yes,noo)
md3 = yy~age + age2 + n.alc
gl3 = glm(md3, family=binomial("logit"))
anova(gl3, test="Chisq")
summary(gl3)

##########################################################################

#Program to analyze the data in Table 11.4 for Example 11.3
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table114.txt")
temp = yd[,1]
oil = yd[,2]
time = yd[,3]
yy = yd[,4]
n = length(yy)
x0= rep(1, n)
#Poisson regression model
md1 = yy~temp + oil + time
gl1 = glm(md1, family=poisson("log"))
anova(gl1, test="Chisq")
summary(gl1)

xx = cbind(x0,temp,oil,time)         
#To find the inital estimates for the regression coefficients
par = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(yy+0.5))

#Fitting Poisson regression model
#f.po is the function that defines the (negative) log likelihood 
f.po <- function(para.p, xx, yy)
{  
m1.x = xx%*%para.p
mu.x = exp(m1.x)
-sum(yy*log(mu.x) - mu.x - lgamma(yy+1))
}
#The function nlminb does not return the hessian matrix
I.p = nlminb(start=par,objective=f.po,xx=xx,yy=yy)
print("Poisson Regression Model")
I.p[[1]] #parameter estimates
I.p[[2]] #minus log-likelihood function
I.p[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.op=nlm(f.po,par,xx=xx,yy=yy,hessian=TRUE)
est=I.op$estimate
log.l=I.op$minimum
hess=I.op$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "temp","oil","time")
res

##########################################################################

#Program to analyze the data in Table 11.4 for NBR model
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table114.txt")
temp = yd[,1]
oil = yd[,2]
time = yd[,3]
yy = yd[,4]
n = length(yy)
x0= rep(1, n)
#Poisson regression model
md1 = yy~temp + oil + time
gl1 = glm(md1, family=poisson("log"))
anova(gl1, test="Chisq")
summary(gl1)

xx = cbind(x0,temp,oil,time)         
#To find the inital estimates for the regression coefficients
par = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(yy+0.5))

#Fitting Poisson regression model
#f.po is the function that defines the (negative) log likelihood 
f.po <- function(para.p, xx, yy)
{  
m1.x = xx%*%para.p
mu.x = exp(m1.x)
-sum(yy*log(mu.x) - mu.x - lgamma(yy+1))
}
#The function nlminb does not return the hessian matrix
I.p = nlminb(start=par,objective=f.po,xx=xx,yy=yy)
print("Poisson Regression Model")
I.p[[1]] #parameter estimates
I.p[[2]] #minus log-likelihood function
I.p[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.op=nlm(f.po,par,xx=xx,yy=yy,hessian=TRUE)
est=I.op$estimate
log.l=I.op$minimum
hess=I.op$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "temp","oil","time")
res

#Fitting negative binomial regression model
#f.nb is the function that defines the (negative) log likelihood 
par.n = rbind(par,0.1)
f.nb <- function(para.n, xx, yy)
{
np=ncol(xx)+1
av=para.n[np]
para.x=para.n[1:(np-1)]  
m1.x = xx%*%para.x
mu.x = exp(m1.x)
de = 1 + av*mu.x
cc = lgamma(yy+1/av) - lgamma(1/av) - lgamma(yy+1)
-sum(cc + yy*log(av*mu.x/de) - log(de)/av)
}
I.nb = nlminb(start=par.n,objective=f.nb,xx=xx,yy=yy)
print("Negative Binomial Regression Model")
I.nb[[1]] #parameter estimates
I.nb[[2]] #minus log-likelihood function
I.nb[[3]] #convergence code: 0 indicates successful

#The function optim returns the hessian matrix with logical TRUE.
I.nb=nlm(f.nb,par.n,xx=xx,yy=yy,hessian=TRUE)
est=I.nb$estimate
log.l=I.nb$minimum
hess=I.nb$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "temp","oil","time","alpha")
res

##########################################################################

#Program to analyze the data in Table 11.4 for GPR model
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\table114.txt")
temp = yd[,1]
oil = yd[,2]
time = yd[,3]
yy = yd[,4]
n = length(yy)
x0= rep(1, n)
#Poisson regression model
md1 = yy~temp + oil + time
gl1 = glm(md1, family=poisson("log"))
anova(gl1, test="Chisq")
summary(gl1)

xx = cbind(x0,temp,oil,time)         
#To find the inital estimates for the regression coefficients
par = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(yy+0.5))

#Fitting Poisson regression model
#f.po is the function that defines the (negative) log likelihood 
f.po <- function(para.p, xx, yy)
{  
m1.x = xx%*%para.p
mu.x = exp(m1.x)
-sum(yy*log(mu.x) - mu.x - lgamma(yy+1))
}
#The function nlminb does not return the hessian matrix
I.p = nlminb(start=par,objective=f.po,xx=xx,yy=yy)
print("Poisson Regression Model")
I.p[[1]] #parameter estimates
I.p[[2]] #minus log-likelihood function
I.p[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.op=nlm(f.po,par,xx=xx,yy=yy,hessian=TRUE)
est=I.op$estimate
log.l=I.op$minimum
hess=I.op$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "temp","oil","time")
res

#Fitting generalized Poisson regression model
#f.gp is the function that defines the (negative) log likelihood 
par.g = rbind(par,0.1)
f.gp <- function(para.g, xx, yy)
{
np=ncol(xx)+1
av=para.g[np]
para.x=para.g[1:(np-1)]  
m1.x = xx%*%para.x
mu.x = exp(m1.x)
ra.x = mu.x/(1+av*mu.x)
de = 1 + av*yy
cc = yy*log(ra.x) + (yy-1)*log(de) - de*ra.x - lgamma(yy+1)
-sum(cc)
}
I.gp = nlminb(start=par.g,objective=f.gp,xx=xx,yy=yy)
print("Generalized Poisson Regression Model")
I.gp[[1]] #parameter estimates
I.gp[[2]] #minus log-likelihood function
I.gp[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.gp=nlm(f.gp,par.g,xx=xx,yy=yy,hessian=TRUE)
est=I.gp$estimate
log.l=I.gp$minimum
hess=I.gp$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "temp","oil","time","alpha")
res

##########################################################################

#Program to analyze domestic violence data in section 11.6
#This prohram estimates ZIP(tau), ZINB(tau) and ZIGP(tau) regression models
library(MASS)
yd=read.table("C:\\Research\\Book\\Data\\section116.txt")
x1 = yd[,1]
x2 = yd[,2]
x3 = yd[,3]
x4 = yd[,4]
x5 = yd[,5]
x6 = yd[,6]
x7 = yd[,7]
x8 = yd[,8]
x9 = yd[,9]
x10= yd[,10]
x11= yd[,11]
x12= yd[,12]
yy = yd[,14]
n = length(yy)
x0= rep(1, n)
#Poisson regression model
md1 = yy~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12
gl1 = glm(md1, family=poisson("log"))
anova(gl1, test="Chisq")
summary(gl1)

xx = cbind(x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)         
#To find the inital estimates for the regression coefficients
par0 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(yy+0.5))
par = rbind(par0, 0.1)

#Fitting zero-inflated Poisson regression model
#fz.po is the function that defines the (negative) log likelihood 
#This log likelihood includes the constant part with no parameters
fz.po <- function(para.p, xx, yy)
{
np=ncol(xx)+1
tau=para.p[np]
para.x=para.p[1:(np-1)]  
m1.x = xx%*%para.x
mu.x = exp(m1.x)
d0=1+mu.x^(-tau)
d1=rep(0,n); d2=rep(0,n);
for (k in 1:n)
{if (yy[k]==0) d1[k]=log(mu.x[k]^(-tau)) + exp(-mu.x[k])
 if (yy[k]>0) d2[k]=yy[k]*log(mu.x[k]) - lgamma(yy[k]+1) - mu.x[k]
}
-sum(-log(d0)+d1+d2)
}
#The function nlminb does not return the hessian matrix
I1.zp = nlminb(start=par,objective=fz.po,xx=xx,yy=yy)
print("Zero Inflated Poisson Regression Model")
I1.zp[[1]] #parameter estimates
I1.zp[[2]] #minus log-likelihood function
I1.zp[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.zp=nlm(fz.po,par,xx=xx,yy=yy,hessian=TRUE)
est=I.zp$estimate
log.l=I.zp$minimum
hess=I.zp$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","tau")
log.l
res
   
#Fitting zero-inflated negative binomial regression model
#fz.nb is the function that defines the (negative) log likelihood 
#This log likelihood includes the constant part with no parameters
fz.nb <- function(para.nb, xx, yy)
{
np=ncol(xx)+2
av=para.nb[np]
tau = para.nb[np-1]
para.x=para.nb[1:(np-2)]  
m1.x = xx%*%para.x
mu.x = exp(m1.x)
d0=1+mu.x^(-tau); dd=av*mu.x; de=1+dd
d1=rep(0,n); d2=rep(0,n); co=rep(0,n)
for (k in 1:n)
{if (yy[k]==0) d1[k]=log(mu.x[k]^(-tau) + de[k]^(-1/av))
 if (yy[k]>0) co[k]=lgamma(yy[k] + 1/av) - lgamma(1/av) - lgamma(yy[k]+1)
 if (yy[k]>0) d2[k]=yy[k]*log(dd[k]/de[k]) - (log(de[k]))/av + co[k]
}
-sum(-log(d0)+d1+d2)
}

#To find the inital estimates for the regression coefficients
par.nb = rbind(par0, 0.1, 0.1)

#The function nlminb does not return the hessian matrix
I1.znb = nlminb(start=par.nb,objective=fz.nb,xx=xx,yy=yy)
print("Zero Inflated Negative Binomial Regression Model")
I1.znb[[1]] #parameter estimates
I1.znb[[2]] #minus log-likelihood function
I1.znb[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
I.znb=nlm(fz.nb,par.nb,xx=xx,yy=yy,hessian=TRUE)
est=I.znb$estimate
log.l=I.znb$minimum
hess=I.znb$hessian
inv=ginv(hess); std=sqrt(diag(inv))
res=cbind(est,std)
colnames(res)=c("estimate", "std. dev.")
rownames(res)=c("constant", "x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","tau","alpha")
log.l
res
  
#Fitting zero-inflated generalized Poisson regression model
#fz.gp is the function that defines the (negative) log likelihood 
#This log likelihood includes the constant part with no parameters
fz.gp <- function(para.gp, xx, yy)
{
np=ncol(xx)+2
av=para.gp[np]
tau = para.gp[np-1]
para.x=para.gp[1:(np-2)]  
m1.x = xx%*%para.x
mu.x = exp(m1.x)
d0=1+mu.x^(-tau); de=1+av*mu.x; dd=1+av*yy
d1=rep(0,n); d2=rep(0,n); co=rep(0,n)
for (k in 1:n)
{if (yy[k]==0) d1[k]=log(mu.x[k]^(-tau) + exp(-mu.x[k]/de[k]))
 if (yy[k]>0) co[k]=lgamma(yy[k]+1) + mu.x[k]*dd[k]/de[k]
 if (yy[k]>0) d2[k]=yy[k]*log(mu.x[k]/de[k]) + (yy[k]-1)*log(dd[k])-co[k]
}
-sum(-log(d0)+d1+d2)
}

#To find the inital estimates for the regression coefficients
par.gp = rbind(par0, 0.1, 0.1)

#The function nlminb does not return the hessian matrix
I1.zgp = nlminb(start=par.gp,objective=fz.gp,xx=xx,yy=yy)
print("Zero Inflated Generalized Poisson Regression Model")
I1.zgp[[1]] #parameter estimates
I1.zgp[[2]] #minus log-likelihood function
I1.zgp[[3]] #convergence code: 0 indicates successful

#The function nlm returns the hessian matrix with logical TRUE.
#This function does not converge for the zero inflated GPR(tau)