I am calculating a likelihood surface for a modified negative binomial distribution with two parameters: mean R and dispersion k (the probability distribution is 'modified' to account for total size distributions of networks, so cant use inherent MLE (i.e. MASS::fitdistr)). Thus, I'm running simulated data through the various functions below to obtain a point estimate and confidence interval.

Since the MLE of k is predicated on R, I am using profile likelihood My approach to this is by creating a matrix of R and k values, then getting the log-likelihood of each permutation and finding where it is maximized.

The current code achieves my goal of calculating a point estimate and confidence interval, but is very time consuming - I don't know if the speed is just inherent with R having to calculate the math for each permutation, or if there is a way to improve the speed through code. With a simulated dataset of length 500, it takes about 30 minutes. I am simulating data from lengths of 50 to 2000 several hundred times times, so imagine some of them will have to run for hours. Improving the speed will save critical days.

The below functions calculate the likelihood surface across a range of R and k values. The function surflike_yn is the time consuming process, as it runs the other functions for each combination of R and k.

My question is how to improve the performance/time of these functions. I typically use R for much simpler processes and have never needed to truly consider performance until now.

#Probability distribution
nb_yn <- function(y,n,R,k) {
  logpyn <- log(n) - log(y) + lgamma(k*y+y-n) - lgamma(k*y) - lgamma(y-n+1) + log((R/k)^(y-n)/(1+R/k)^(k*y+y-n))

# Likelihood function
like_yn<- function(simdata,r1,k1){
  runs <- nrow(simdata)
  liks <- rep(1,runs)
  for(z in 1:runs){
    liks[z] <- nb_yn(simdata[z,1],simdata[z,2],r1,k1)

# Calculate likelihood surface (***this is the time consuming step***)
surflike_yn<- function(simdata,Rrange,krange){
  likesurf <- data.frame(matrix(NA, nrow=length(Rrange),length(krange)))
  for(i in 1:length(Rrange)){
    for(j in 1:length(krange)){
      likesurf[i,j] <- like_yn(simdata,Rrange[i],krange[j])

#Get point and CI estimates
  chiV<-qchisq(conf.interval/100, df=1)/2 
  prfk <- apply(likesurf,2,function(x){max(x)})
  prfk2 <- krange[prfk-max(prfk)>-chiV]
  prfR <- apply(likesurf,1,function(x){max(x)})
  prfR2 <- Rrange[prfR-max(prfR)>-chiV]

  output <- matrix(NA,2,3)
  output[1,1] <- Rrange[sum(seq(1,length(Rrange))%*%likesurf_max)]
  output[1,2] <- min(prfR2)
  output[1,3] <- max(prfR2)
  output[2,1] <- krange[sum(likesurf_max%*%seq(1,length(krange)))]
  output[2,2] <- min(prfk2)
  output[2,3] <- ifelse(max(prfk2)==max(krange),Inf,max(prfk2))
  colnames(output) <- c("point_est","lower_ci","upper_ci"); rownames(output) <- c("R","k")

Below, I am providing code that I use to simulate data to run these functions (this does not necessarily need review; I use these functions in other processes which is why they may not be optimized)

My simulated data are a two column matrix. The first column are the size distributions (Y), the second is how many values this is conditioned on (i.e. P(Y|n) if two networks cant be unambiguously teases apart - n is almost always 1 in practice)

#Set up R and k "search grid" for profile likelihood
Rlow=0.05; Rhigh=1.05; Rrange=seq(Rlow,Rhigh,by=0.01)
klow=0.04; khigh=55; krange=seq(klow,khigh,by=0.01)

#Generate dummy data with R=0.5 and k=0.25 and ALL index cases=1
  #Brancing process function
bp <- function(gens=20, init.size=1, offspring, ...){  
  Z <- list(); Z[[1]] <- init.size; i <- 1 
  while(sum(Z[[i]]) > 0 && i <= gens) { 
    Z[[i+1]] <- offspring(sum(Z[[i]]), ...) 
    i <- i+1 
  #generate data (n=100)
ind<-replicate(100,bp(offspring=rnbinom, size=0.25,mu=0.5)) 
clust<-t(data.frame(lapply(ind,function(x) sum(unlist(x)))))
Yn <- cbind(clust,rep(1,times=length(clust))) #2 x n matrix will all index cases=1

#Calculate surface likelihoods and parameter estimates
  start_time <- Sys.time()
surf <- surflike_yn(Yn,Rrange,krange)
  end_time <- Sys.time()
  duration <- end_time - start_time
surf_max <- surf==max(surf)
  • \$\begingroup\$ Hi @greybeard - like_yn() should be the second function listed in the first code block. Please let me know if it isn't showing up, it appears on my end. Thank you! \$\endgroup\$
    – jpsmith
    May 5, 2020 at 11:56
  • 1
    \$\begingroup\$ (By the looks of it, I didn't copy from code block to my IDE quite successfully.) \$\endgroup\$
    – greybeard
    May 5, 2020 at 20:00


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.