Как рассчитать BIC и AIC для модели gmm в R с использованием plm?

Я оцениваю модель GMM, используя библиотеку plm. У меня разные моментные условия.

Z <- list(~YDWPP + ST_DEGREE, ~YDWPP + ST_DEGREE, ~YDWPP + ST_DEGREE, 
    ~YDWPP + ST_DEGREE, ~YDWPP + ST_TRANSITIVITY, ~YDWPP + ST_STRUC_HOLE, 
    ~YDWPP + ST_STRUC_HOLE, ~YDWPP + ST_STRUC_HOLE, ~YDWPP + 
        ST_STRUC_HOLE)

Z <- lapply(Z, as.formula)

lg.gmm <- list(c(4L, 8L), c(5L, 8L), c(6L, 8L), 7:8, 7:8, c(4L, 8L), c(5L, 
8L), c(6L, 8L), 7:8)

Я запускаю цикл для каждого набора ограничений момента Z, так что

out.1 <- list()
for(i in seq_along(Z)){
  plm.gmm <-
  pgmm(
  dynformula(as.formula(model), lg),
  data = pdata,
  effect = 'twoway',
  model = 'twostep',
  transformation = 'd',
  gmm.inst = Z[[i]],
  lag.gmm =  c(lg.gmm[[i]][[1]], lg.gmm[[i]][[2]])
  )
sum <- summary(plm.gmm, robust = T)
print(sum)
out.1[[i]] <- sum
}

Я хотел бы сравнить эти модели, используя, например, BIC и AIC

AIC(plm.gmm, k=2)
Error in UseMethod("logLik") : 
  no applicable method for 'logLik' applied to an object of class "c('pgmm', 'panelmodel')"

Любые идеи о том, как вычислить BIC и AIC или альтернативные методы выбора между различными ограничениями момента?


r benchmarking plm gmm
person Mario GS    schedule 12.09.2017    source источник
comment
Привет, нашел решение проблемы? Я застрял на той же ошибке.   -  person Daniel    schedule 07.12.2017

Ответы (2)


arrow_upward
4
arrow_downward

Я слежу за ответом на этот вопрос.

Дополнительные сведения о критериях AIC можно найти в Википедии.

Вот код, который должен работать. Однако вы не предоставили никакой воспроизводимой оценки модели. Следовательно, это без проверки для вашего случая.

# Function: Calculates AIC based on an lm or plm object

AIC_adj <- function(mod){
  # Number of observations
  n.N   <- nrow(mod$model)
  # Residuals vector
  u.hat <- residuals(mod)
  # Variance estimation
  s.sq  <- log( (sum(u.hat^2)/(n.N)))
  # Number of parameters (incl. constant) + one additional for variance estimation
  p     <-  length(coef(mod)) + 1

  # Note: minus sign cancels in log likelihood
  aic <- 2*p  +  n.N * (  log(2*pi) + s.sq  + 1 ) 

  return(aic)
}
person Alex    schedule 13.04.2018
comment
Может быть, стоит отметить одну вещь? Я просто застрял из-за этого: эффект = «двусторонняя» модель plm, и формула Алекса не будет включать время и отдельные эффекты в «p» (количество параметров) здесь. В вашем исходном вопросе вы могли бы написать фиктивную регрессию, а затем AIC() включил бы эти манекены в «p». Обычно вы, вероятно, не хотите этого, но все же важно убедиться, что мы сравниваем. - person Jakob; 03.09.2019

arrow_upward
1
arrow_downward

Необходимо учитывать разные габариты (и количество параметров) различных вариантов моделей панелей. Продолжая предыдущий пример:

    aicbic_plm <- function(object, criterion) {


    # object is "plm", "panelmodel" 
    # Lets panel data has index :index = c("Country", "Time")

    sp = summary(object)

    if(class(object)[1]=="plm"){
    u.hat <- residuals(sp) # extract residuals
    df <- cbind(as.vector(u.hat), attr(u.hat, "index"))
    names(df) <- c("resid", "Country", "Time")
    c = length(levels(df$Country)) # extract country dimension 
    t = length(levels(df$Time)) # extract time dimension 
    np = length(sp$coefficients[,1]) # number of parameters
    n.N = nrow(sp$model) # number of data
    s.sq  <- log( (sum(u.hat^2)/(n.N))) # log sum of squares

    # effect = c("individual", "time", "twoways", "nested"),
    # model = c("within", "random", "ht", "between", "pooling", "fd")

    # I am making example only with some of the versions:

    if (sp$args$model == "within" & sp$args$effect == "individual"){
    n = c
    np = np+n+1 # update number of parameters
    }

    if (sp$args$model == "within" & sp$args$effect == "time"){
    T = t
    np = np+T+1 # update number of parameters
    }

    if (sp$args$model == "within" & sp$args$effect == "twoways"){
    n = c
    T = t
    np = np+n+T # update number of parameters
    }
    aic <- round(       2*np  +  n.N * (  log(2*pi) + s.sq  + 1 ),1)
    bic <- round(log(n.N)*np  +  n.N * (  log(2*pi) + s.sq  + 1 ),1)

    if(criterion=="AIC"){
    names(aic) = "AIC"
    return(aic)
    }
    if(criterion=="BIC"){
    names(bic) = "BIC"
    return(bic)
    }
    }
    }
person Rookie    schedule 07.12.2018