[R] Fitting a Bivariate Poisson log-normal distribution

From: Famoye, Felix <famoy1kf_at_cmich.edu>
Date: Sun, 27 Jul 2008 15:36:41 -0400


I wrote an R program to fit a bivariate Poisson log-normal model. The model is first proposed by Aitchison and Ho (1989) in Biometrika, Volume 76 pages 643-653. The model involves using a double integration. The following is the program and the data that I used. My major problem was that R was running for a long time (more than 3 hours) with no results. I would like to place some printing at some places to see if R was doing something. Unfortunately, there was nothing printed when I included printing in the function "f". Note that I have used the package "adapt" for double integration. Can someone provide some help on this problem? I would like to know if R is doing anything at all in fitting the model. Thank you for your time.

#This program will be used to estimate the parameters of
Poisson-lognormal model.
yd = read.table(file="H:\\Bivariate\\bact.txt")  

y1 = yd[[1]]
y2 = yd[[2]]
n = length(y1)
x0= rep(1, n)

y1  = cbind(y1)
y2  = cbind(y2)
yy=cbind(y1,y2)
xx = cbind(x0,yd[[3]])    
av1=mean(y1)

av2=mean(y2)
cov=var(yy)
mu=cbind(av1, av2)
rho=cov[1,2]/sqrt(cov[1,1]*cov[2,2])
sd1=sqrt(cov[1,1])
sd2=sqrt(cov[2,2]) 
oth=rbind(sd1,sd2,rho)

z = vector(length=2) #This is a column vector with two rows

#To find the inital estimates for the regression coefficients

para1 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(y1+0.5))
para2 = ginv(t(xx) %*% xx) %*% (t(xx) %*% log(y2+0.5))
parab = rbind(para1, para2, oth)

parab

#---------------------------------------------------------------------
# Fitting Poisson lognorma model
#f is the function that defines the (negative) log likelihood
f <- function(parab,y1,y2,xx)
{
vv1=(parab[5])^2
vv2=(parab[6])^2
para.1= cbind(parab[1:2])
para.2= cbind(parab[3:4])
mu.1 = xx%*%para.1 - vv1/2
#mu.1 = exp(mu.1)
mu.2 = xx%*%para.2 - vv2/2
#mu.2 = exp(mu.2)

mu.v=cbind(mu.1,mu.2)

cov[1,1]=vv1
cov[2,2]=vv2
cov[1,2]=parab[7]*parab[5]*parab[6]
cov[2,1]=cov[1,2]

n=length(y1)
va=rep(0.0, n)
for (i in 1:n)
{
fred=function(z)
{
tz=log(z)-mu.v[i,]
qq=(t(tz)%*%ginv(cov))%*%tz
fq=(exp(-qq/2))/(2*pi*z[1]*z[2]*sqrt(det(cov)))
p1=((z[1]^y1[i])*exp(-z[1]))/gamma(y1[i]+1)
p2=((z[2]^y2[i])*exp(-z[2]))/gamma(y2[i]+1)
int=fq*p1*p2
}
va[i]=adapt(2,lo=c(0,0),up=c(100,100),functn=fred,min=1000,eps=1.e-6,mu. v=mu.v,cov=cov,y1=y1,y2=y2)$value
}
-sum(va)
}

I.p = nlminb(start=parab,objective=f,y1=y1,y2=y2,xx=xx) print("Poisson Log-normal Regression Model") I.p

   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   0 0 0
   1 0 0
   1 0 0
   1 0 0
   1 0 0
   1 0 0
   1 0 0
   1 0 0
   2 0 0
   2 0 0
   2 0 0
   2 0 0
   2 0 0
   2 0 0
   3 0 0
   3 0 0
   3 0 0
   3 0 0
   3 0 0
   3 0 0
   4 0 0
   4 0 0
   4 0 0
   5 0 0
   5 0 0
   7 0 0
   8 0 0
   9 0 0
   9 0 0
  12 0 0
  13 0 0
  14 0 0
  16 0 0
  20 0 0
   0 1 0
   0 1 0
   0 1 0
   0 1 0
   0 1 0
   0 1 0
   0 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   1 1 0
   2 1 0
   2 1 0
   2 1 0
   2 1 0
   2 1 0
   2 1 0
   2 1 0
   3 1 0
   3 1 0
   3 1 0
   4 1 0
   4 1 0
   5 1 0
   5 1 0
   5 1 0
   9 1 0
  10 1 0
  10 1 0
  15 1 0
   0 2 0
   0 2 0
   0 2 0
   0 2 0
   0 2 0
   0 2 0
   0 2 0
   0 2 0
   0 2 0
   1 2 0
   1 2 0
   1 2 0
   1 2 0
   1 2 0
   1 2 0
   1 2 0
   1 2 0
   1 2 0
   2 2 0
   2 2 0
   2 2 0
   3 2 0
   4 2 0
   4 2 0
   5 2 0
   6 2 0
   6 2 0
   6 2 0
   7 2 0
   7 2 0
  12 2 0
  15 2 0
  26 2 0
   0 3 0
   0 3 0
   0 3 0
   0 3 0
   0 3 0
   0 3 0
   0 3 0
   0 3 0
   0 3 0
   1 3 0
   1 3 0
   1 3 0
   1 3 0
   1 3 0
   2 3 0
   2 3 0
   4 3 0
   5 3 0
   7 3 0
  10 3 0
   0 4 0
   0 4 0
   0 4 0
   0 4 0
   1 4 0
   1 4 0
   1 4 0
   1 4 0
   3 4 0
   3 4 0
   0 5 0
   0 6 0
   0 6 0
   1 6 0
   3 6 0
   0 7 0
   0 7 0
   1 7 0
   0 8 0
   1 8 0
   0 9 0
   0 9 0
   0 9 0
   0 0 1
   0 0 1
   0 0 1
   0 0 1
   0 0 1
   0 0 1
   0 0 1
   0 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   1 0 1
   2 0 1
   2 0 1
   2 0 1
   2 0 1
   2 0 1
   2 0 1
   3 0 1
   3 0 1
   3 0 1
   3 0 1
   3 0 1
   3 0 1
   3 0 1
   3 0 1
   4 0 1
   4 0 1
   4 0 1
   4 0 1
   4 0 1
   4 0 1
   4 0 1
   4 0 1
   5 0 1
   5 0 1
   7 0 1
   7 0 1
   8 0 1
   9 0 1
   9 0 1
  11 0 1
  15 0 1
  15 0 1
  19 0 1
  20 0 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   0 1 1
   1 1 1
   1 1 1
   1 1 1
   2 1 1
   2 1 1
   2 1 1
   2 1 1
   2 1 1
   4 1 1
   4 1 1
   5 1 1
   5 1 1
   5 1 1
   6 1 1
   8 1 1
   8 1 1
   8 1 1
   8 1 1
   9 1 1
  23 1 1
  23 1 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   0 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   1 2 1
   2 2 1
   2 2 1
   2 2 1
   3 2 1
   3 2 1
   3 2 1
   4 2 1
   4 2 1
   6 2 1
   0 3 1
   0 3 1
   0 3 1
   0 3 1
   0 3 1
   0 3 1
   0 3 1
   1 3 1
   1 3 1
   1 3 1
   1 3 1
   1 3 1
   2 3 1
   2 3 1
   2 3 1
   3 3 1
   6 3 1
   0 4 1
   0 4 1
   0 4 1
   0 4 1
   0 4 1
   1 4 1
   1 4 1
   1 4 1
   1 4 1
   2 4 1
   5 4 1
   0 5 1
   0 5 1
   0 5 1
   1 5 1
   1 5 1
   0 6 1
   0 6 1
   1 7 1
   1 7 1
   1 8 1
   0 10 1
   0 10 1




Felix Famoye
Department of Mathematics
Central Michigan University
Mt. Pleasant, Michigan 48859-0001

e-mail: felix.famoye_at_cmich.edu
web site: http://www.cst.cmich.edu/users/famoy1kf/ voice: (989)774-5497, fax: (989)774-2414



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 27 Jul 2008 - 20:06:49 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 27 Jul 2008 - 21:33:08 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive