[R] Fitting a Bivariate Poisson log-normal distribution

Famoye, Felix famoy1kf at cmich.edu
Sun Jul 27 21:36:41 CEST 2008


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



More information about the R-help mailing list