Средние предельные эффекты (AME) в модели с частичными пропорциональными шансами

Question

Как получить средние предельные эффекты (AME) для каждой категории / порога в модели с частичным пропорциональным коэффициентом (PPOM)?

Это мой первый пост на этом форуме. Я надеюсь, что прислушался к самым важным рекомендациям, чтобы задавать хорошие вопросы.

Этот выборочный набор данных состоит из порядковой переменной результата (Y1) и трех независимых переменных (VAR1, VAR2, VAR3).

set.seed(3)
sampleData <- data.frame(id = 1:1000, Y1 = sample(c("1", "2", "3", "4"), 
    1000, replace=TRUE), Var1 = rnorm(1000, 40, 10), 
    Var2 = rnorm(1000, 60, 10), Var3 = rnorm(1000, 80, 5))

Предполагая, что допущение пропорциональных коэффициентов нарушено, можно выполнить модель частичных пропорциональных коэффициентов (PPOM), используя пакет ordinal для прогнозирования Y1 тремя независимыми переменными (Var1, Var2, Var3).

library(ordinal)  
PPOM <- clm(as.factor(Y1) ~ Var1 + Var2 + Var3, 
        nominal = ~ Var1 + Var2 + Var3, data = sampleData)

Мы получаем следующий вывод с коэффициентами для каждой категории:

summary(PPOM)

formula: as.factor(Y1) ~ Var1 + Var2 + Var3
nominal: ~Var1 + Var2 + Var3
data:    sampleData

 link  threshold nobs logLik   AIC     niter max.grad cond.H 
 logit flexible  1000 -1381.17 2786.34 4(0)  2.82e-10 2.2e+07

Coefficients: (3 not defined because of singularities)
     Estimate Std. Error z value Pr(>|z|)
Var1       NA         NA      NA       NA
Var2       NA         NA      NA       NA
Var3       NA         NA      NA       NA

Threshold coefficients:
                  Estimate Std. Error z value
1|2.(Intercept)  0.4952642  1.2260010   0.404
2|3.(Intercept)  1.9790234  1.0724982   1.845
3|4.(Intercept)  2.0892425  1.2550636   1.665
1|2.Var1         0.0026194  0.0075920   0.345
2|3.Var1        -0.0077578  0.0065845  -1.178
3|4.Var1        -0.0064243  0.0075364  -0.852
1|2.Var2        -0.0001089  0.0074568  -0.015
2|3.Var2        -0.0082836  0.0063447  -1.306
3|4.Var2        -0.0073638  0.0071008  -1.037
1|2.Var3        -0.0219767  0.0140701  -1.562
2|3.Var3        -0.0157235  0.0121943  -1.289
3|4.Var3        -0.0047098  0.0141844  -0.332

Меня интересуют AME для каждого предиктора для каждой категории. Используя margins, я получаю только AME для всех порогов.

library(margins) 
summary(margins(PPOM))

Выход:

 factor    AME     SE      z      p   lower  upper
   Var1 0.0000 0.0000 1.1365 0.2557 -0.0000 0.0001
   Var2 0.0000 0.0000 1.3056 0.1917 -0.0000 0.0001
   Var3 0.0001 0.0001 0.9990 0.3178 -0.0001 0.0002

Кто-нибудь знает, как рассчитать AME для каждой категории?

Любая помощь будет принята с благодарностью!

Answer 1

Я знаю, что немного опоздал, но я только что написал короткую (и довольно неэффективную) программу для ручного расчета средних предельных эффектов для этого типа модели.

Средние предельные эффекты для частично-пропорционального порядкового логита / пробита рассчитываются так же, как и для обычного порядкового логита / пробита. Ниже приведен отрывок из «Микроэконометрики Камерона и Триведи: методы и приложения». Предельные эффекты для непрерывных переменных задаются уравнением в нижней части экрана. Предельный эффект для дискретной переменной равен Pr (y_i = j | x = 1) - Pr (y_i = j | x = 0); Pr (y_i = j) задается в уравнении в верхней части страницы.

Я бы предложил вручную рассчитать предельные эффекты для каждого набора оцениваемых бета-версий. Это включает в себя сохранение матриц / векторов, которые состоят из ваших данных, ваших разных бета-версий и коэффициента в переменной, для которой вы рассчитываете ME. Это также означает, что вы вычисляете разные ME для каждого результата и каждого коэффициента (в вашем случае 3 * число порядковых результатов).

Ниже приведен код, который я написал, чтобы сделать это сам. Это не эффективно и идет один за другим. Он рассчитывает значения ME фиктивной переменной под названием «snap12» для трех результатов, поэтому существует шесть уравнений. «Cut1» и «Cut2» - это точки среза (обычно называемые альфа в уравнении ME). Xdat - это данные, исключая переменную, для которой я вычисляю ME. xcoef1 - это первый набор коэффициентов (за исключением переменной ME), а xcoef2 - это второй набор коэффициентов. Snapcoef1 и snapcoef2 - это коэффициенты для переменной ME. Наконец, pnorm дает значение CDF, которое требуется для порядкового пробита ME.

Надеюсь, это кому-нибудь поможет!

##for pr(v low secure)

#1st set of coefficients
p3a = ((1 - pnorm(cut2 - xdat%*%xcoef1 - snapcoef1)) - 
       (1 - pnorm(cut2 - xdat%*%xcoef1)))
#2nd set
p3b = ((1 - pnorm(cut2 - xdat%*%xcoef2 - snapcoef2)) - 
       (1 - pnorm(cut2 - xdat%*%xcoef2)))

#for pr(low sec)
#1st set
p2a = (pnorm(cut2 - xdat%*%xcoef1 - snapcoef1) - pnorm(cut1 - xdat%*%xcoef1 - snapcoef1)) - 
         (pnorm(cut2 - xdat%*%xcoef1) - pnorm(cut1 - xdat%*%xcoef1))

p2b = (pnorm(cut2 - xdat%*%xcoef2 - snapcoef2) - pnorm(cut1 - xdat%*%xcoef2 - snapcoef2)) - 
  (pnorm(cut2 - xdat%*%xcoef2) - pnorm(cut1 - xdat%*%xcoef2))

## for pr(fodo sec)

p1a = (pnorm(cut1 - xdat%*%xcoef1 - snapcoef1)) - (pnorm(cut1 - xdat%*%xcoef1))

p1b = (pnorm(cut1 - xdat%*%xcoef2 - snapcoef2)) - (pnorm(cut1 - xdat%*%xcoef2))