Log Likelihood Calculations for Censored Data in R
Your function p_function is vectorized and that is nice. There's a little difficulty for the sum though.
Below is the improved code that takes advantage of the vectorization. I also replaced log(1+x) with log1p(x), which is better for a small |x|. And I replaced the stuff with the lgamma functions with some stuff with the lbeta function.
new_likelihood2 &lt;- function(Y, R, k) {

  p_function &lt;- function(y, n){ 
    exp(log(n) - log(y) - lbeta(k*y, y-n+1) - log(k*y+y-n) + (y-n)*log(R/k) - (k*y+y-n)*log1p(R/k))
  }
  
  p_function_sum &lt;- function(y, n){
      vapply(
        seq_along(y), function(i) sum(p_function(1:(y[i]-1), n[i])), numeric(1)
      )
  }
  
  uncensored &lt;- Y[, 3] == 0
  Yuncensored &lt;- Y[uncensored, 1:2]
  Ycensored   &lt;- Y[!uncensored, 1:2]
  liks_a &lt;- log(p_function(Yuncensored[, 1], Yuncensored[, 2]))
  liks_b &lt;- log1p(-p_function_sum(Ycensored[, 1], Ycensored[, 2]))
  
  sum(liks_a, liks_b)
}