Error in as.vector(x, mode) : cannot coerce type 'closure' to vector of type 'any' -- when running a nested function
get_pi_ij <- function() {
  alphas <-
    matrix(get_alpha(),
           nrow = length(get_alpha()),
           ncol = length(get_alpha()),
           byrow = FALSE)
  betas <-
    matrix(get_beta_prelim(), # your were missing "()"
           nrow = length(get_beta_prelim()),
           ncol = length(get_beta_prelim()),
           byrow = TRUE) 
  pi_ij <- exp(alphas + betas + gamma * c)
  return(pi_ij)
}


Problem: running get_pi_ij() gives the error:
Error in as.vector(x, mode) : cannot coerce type 'closure' to vector of type 'any'
Called from: as.vector(data)
The first thing this function does is to make the resulting alphas and beta_prelims into matrixes that match c so that they can be calculated together. This is where something goes wrong, and I have not been able to figure out what. If I use <<- to save the alphas and betas to the global environment in the prior functions for alphas and betas and replace that with those in the faulty function, it works. So I assume it has to do with how I call the functions inside the matrix creation.
get_pi_ij <- function() {
  alphas <-
    matrix(get_alpha(),
           nrow = length(get_alpha()),
           ncol = length(get_alpha()),
           byrow = FALSE)
  betas <-
    matrix(get_beta_prelim,
           nrow = length(get_beta_prelim()),
           ncol = length(get_beta_prelim()),
           byrow = TRUE) 
  pi_ij <- exp(alphas + betas + gamma * c)
  return(pi_ij)
}
get_pi_ij()

I added the full code cause it's not too long and the first parts are just definitions. Makes it easier to test it.
Everything up to the final function works as it is supposed to
size <- 18 
gamma <- -0.07 

c <- structure(c(0, 4, 8, 12, 16, 4, 4, 8, 12, 16, 8, 8, 8, 12, 16, 
            12, 12, 16, 16, 4, 0, 4, 8, 12, 8, 4, 4, 8, 12, 12, 8, 8, 12, 
            16, 12, 12, 16, 16, 8, 4, 0, 4, 8, 12, 8, 4, 8, 12, 16, 12, 8, 
            12, 16, 16, 12, 16, 16, 12, 8, 4, 0, 4, 12, 8, 4, 4, 8, 16, 12, 
            8, 8, 12, 16, 12, 12, 16, 16, 12, 8, 4, 0, 16, 12, 8, 4, 4, 16, 
            12, 8, 8, 8, 16, 12, 12, 16, 4, 8, 12, 12, 16, 0, 4, 8, 12, 16, 
            4, 4, 8, 12, 16, 8, 8, 12, 12, 4, 4, 8, 8, 12, 4, 0, 4, 8, 12, 
            8, 4, 4, 8, 12, 8, 8, 12, 12, 8, 4, 4, 4, 8, 8, 4, 0, 4, 8, 12, 
            8, 4, 8, 12, 12, 8, 12, 12, 12, 8, 8, 4, 4, 12, 8, 4, 0, 4, 12, 
            8, 4, 4, 8, 12, 8, 8, 12, 16, 12, 12, 8, 4, 16, 12, 8, 4, 0, 
            16, 12, 8, 4, 4, 12, 8, 8, 12, 8, 12, 16, 16, 16, 4, 8, 12, 12, 
            16, 0, 4, 8, 12, 16, 4, 8, 12, 8, 8, 8, 12, 12, 12, 4, 4, 8, 
            8, 12, 4, 0, 4, 8, 12, 4, 4, 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 
            8, 8, 4, 0, 4, 8, 8, 4, 8, 8, 12, 12, 12, 8, 8, 12, 8, 8, 4, 
            4, 12, 8, 4, 0, 4, 8, 4, 4, 8, 16, 16, 16, 12, 8, 16, 12, 12, 
            8, 4, 16, 12, 8, 4, 0, 12, 8, 4, 8, 12, 12, 16, 16, 16, 8, 8, 
            12, 12, 12, 4, 4, 8, 8, 12, 0, 4, 8, 4, 12, 12, 12, 12, 12, 8, 
            8, 8, 8, 8, 8, 4, 4, 4, 8, 4, 0, 4, 4, 16, 16, 16, 12, 12, 12, 
            12, 12, 8, 8, 12, 8, 8, 4, 4, 8, 4, 0, 4, 16, 16, 16, 16, 16, 
            12, 12, 12, 12, 12, 8, 8, 8, 8, 8, 4, 4, 4, 0), 
          .Dim = c(19L, 19L), 
          .Dimnames = list(NULL, c("X1", "X2", "X3", "X4", "X5", 
                                   "X6", "X7", "X8", "X9", "X10", 
                                   "X11", "X12", "X13", "X14", "X15", 
                                   "X16", "X17", "X18", "X19")))

set.seed(12345) 
h_share <- diff(c(0, sort(runif(size)), 1)) 
e_share <- diff(c(0, sort(runif(size)), 1)) 

alpha <- numeric(size + 1)
beta <- numeric(size + 1)

get_beta_prelim <- function() {
  a_matrix <- exp(alpha + t(c) * gamma)
  beta_prelim <- log(e_share) - log(colSums(a_matrix))
  return(beta_prelim)
}

get_beta <- function() {
  beta <- get_beta_prelim() - get_beta_prelim()[[1]]
  return(beta)
}

get_alpha_prelim <- function() {
  b_matrix <- t(exp(beta + t(c) * gamma))
  alpha_prelim <- log(h_share) - log(rowSums(b_matrix))
  return(alpha_prelim)
}

get_alpha <- function() {
  alpha <- get_alpha_prelim() - get_alpha_prelim()[[1]]
  return(alpha)
}

get_pi_ij <- function() {
  alphas <-
    matrix(get_alpha(),
           nrow = length(get_alpha()),
           ncol = length(get_alpha()),
           byrow = FALSE)
  betas <-
    matrix(get_beta_prelim,
           nrow = length(get_beta_prelim()),
           ncol = length(get_beta_prelim()),
           byrow = TRUE) 
  pi_ij <- exp(alphas + betas + gamma * c)
  return(pi_ij)
}
get_pi_ij()

