ช่วงความเชื่อมั่นสำหรับโมเดลผสมที่มีอุปสรรคหรือสูงเกินจริงเป็นศูนย์

ฉันต้องการคำนวณ CI ในแบบจำลองแบบผสม แบบจำลองทวินามลบที่พองตัวเป็นศูนย์ และแบบจำลองอุปสรรค์ รหัสของฉันสำหรับโมเดลอุปสรรค์มีลักษณะดังนี้ (x1, x2 ต่อเนื่อง, หมวดหมู่ x3):

m1 <- glmmTMB(count~x1+x2+x3+(1|year/class),
          data = bd, zi = ~x2+x3+(1|year/class), family = truncated_nbinom2,
          )

ฉันใช้ confint และได้ผลลัพธ์เหล่านี้:

ci <- confint(m1,parm="beta_")
ci
                                          2.5 %       97.5 %     Estimate
cond.(Intercept)                   1.816255e-01  0.448860094  0.285524861
cond.x1                            9.045278e-01  0.972083366  0.937697401
cond.x2                            1.505770e+01 26.817439186 20.094998772
cond.x3high                        1.190972e+00  1.492335046  1.333164894
cond.x3low                         1.028147e+00  1.215828654  1.118056377
cond.x3reg                         1.135515e+00  1.385833853  1.254445909
class:year.cond.Std.Dev.(Intercept)2.256324e+00  2.662976154  2.441845815
year.cond.Std.Dev.(Intercept)      1.051889e+00  1.523719169  1.157153015
zi.(Intercept)                     1.234418e-04  0.001309705  0.000402085
zi.x2                              2.868578e-02  0.166378014  0.069084606
zi.x3high                          8.972025e-01  1.805832900  1.272869874

ฉันคำนวณช่วงเวลาได้อย่างถูกต้องหรือไม่? เหตุใดจึงมีเพียงหมวดหมู่เดียวใน x3 สำหรับ zi หากเป็นไปได้ ฉันอยากจะทราบว่าเป็นไปได้หรือไม่ที่จะพล็อต CI เหล่านี้

ขอบคุณ!

ข้อมูลมีลักษณะดังนี้:

     class id year count  x1      x2        x3 
956    5 3002 2002      3 15.6    47.9      high
957    5 4004 2002      3 14.3    47.9      low
958    5 6021 2002      3 14.2    47.9      high
959    4 2030 2002      3 10.5    46.3      high
960    4 2031 2002      3 15.3    46.3      high
961    4 2034 2002      3 15.2    46.3      reg

ด้วย x1 และ x2 ต่อเนื่อง, x3 ตัวแปรหมวดหมู่สามระดับ (ตัวประกอบ)

สรุปโมเดล:

summary(m1)
'giveCsparse' has been deprecated; setting 'repr = "T"' for you'giveCsparse' has been deprecated; setting 'repr = "T"' for you'giveCsparse' has been deprecated; setting 'repr = "T"' for you
 Family: truncated_nbinom2  ( log )
Formula:          count ~ x1 + x2 + x3 + (1 | year/class)
Zero inflation:          ~x2 + x3 + (1 | year/class)
Data: bd

     AIC      BIC   logLik deviance df.resid 
 37359.7  37479.7 -18663.8  37327.7    13323 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 class:year(Intercept) 0.79701  0.8928  
 year      (Intercept) 0.02131  0.1460  
Number of obs: 13339, groups:  class:year, 345; year, 15

Zero-inflation model:
 Groups    Name        Variance  Std.Dev. 
 dpto:year (Intercept) 1.024e+02 1.012e+01
 year      (Intercept) 7.842e-07 8.856e-04
Number of obs: 13339, groups:  class:year, 345; year, 15

Overdispersion parameter for truncated_nbinom2 family (): 1.02 

Conditional model:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -1.25343    0.23081  -5.431 5.62e-08 ***
x1             -0.06433    0.01837  -3.501 0.000464 ***
x2              3.00047    0.14724  20.378  < 2e-16 ***
x3high          0.28756    0.05755   4.997 5.82e-07 ***
x3low           0.11159    0.04277   2.609 0.009083 ** 
x3reg           0.22669    0.05082   4.461 8.17e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -7.8188     0.6025 -12.977  < 2e-16 ***
x2              -2.6724     0.4484  -5.959 2.53e-09 ***
x3high           0.2413     0.1784   1.352  0.17635    
x3low           -0.1325     0.1134  -1.169  0.24258    
x3reg           -0.3806     0.1436  -2.651  0.00802 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

CI กับไม้กวาดผสม

> broom.mixed::tidy(m1, effects="fixed", conf.int=TRUE)
# A tibble: 12 x 9
   effect component term           estimate std.error statistic  p.value conf.low conf.high
   <chr>  <chr>     <chr>             <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
 1 fixed  cond      (Intercept)     -1.25      0.231      -5.43 5.62e- 8  -1.71     -0.801 
 2 fixed  cond      x1              -0.0643    0.0184     -3.50 4.64e- 4  -0.100    -0.0283
 3 fixed  cond      x2               3.00      0.147      20.4  2.60e-92   2.71      3.29  
 4 fixed  cond      x3high           0.288     0.0575      5.00 5.82e- 7   0.175     0.400 
 5 fixed  cond      x3low            0.112     0.0428      2.61 9.08e- 3   0.0278    0.195 
 6 fixed  cond      x3reg            0.227     0.0508      4.46 8.17e- 6   0.127     0.326 
 7 fixed  zi        (Intercept)     -9.88      1.32       -7.49 7.04e-14 -12.5      -7.30  
 8 fixed  zi        x1               0.214     0.120       1.79 7.38e- 2  -0.0206    0.448 
 9 fixed  zi        x2              -2.69      0.449      -6.00 2.01e- 9  -3.57     -1.81  
10 fixed  zi        x3high           0.232     0.178       1.30 1.93e- 1  -0.117     0.582 
11 fixed  zi        x3low           -0.135     0.113      -1.19 2.36e- 1  -0.357     0.0878
12 fixed  zi        x4reg           -0.382     0.144      -2.66 7.74e- 3  -0.664    -0.101 

person MKie45    schedule 12.04.2021    source แหล่งที่มา
comment
(1) ใช่ เท่าที่ฉันสามารถบอกได้ (2) ยากที่จะพูดหากไม่มี ตัวอย่างที่สามารถทำซ้ำได้น้อยที่สุด ฉันยอมรับว่ามันดูแปลก (3) คุณสามารถใช้ broom.mixed::tidy(m1, effects="fixed", conf.int=TRUE) เพื่อรับข้อมูลในรูปแบบที่สะดวกสำหรับ ggplot2 (หรือแม้แต่ dotwhisker::dwplot(m1, effects="fixed"))   -  person Ben Bolker    schedule 12.04.2021
comment
ฉันเพิ่มตัวอย่างข้อมูลและคำอธิบายของตัวแปร ปัญหาคือ x3 มีสามหมวดหมู่ แต่ confint รายงานเพียงหมวดหมู่เดียว ยังไงก็จะลองใช้ broom.mixed ดูครับ   -  person MKie45    schedule 12.04.2021
comment
ฉันจะลองดูว่าฉันจะไปได้ไกลกว่านี้หรือเปล่า แต่ฉันก็ไม่แน่ใจ ฉัน/เราต้องการตัวอย่างที่ ทำซ้ำได้ (ดูลิงก์ ตัวอย่างที่ทำซ้ำได้น้อยที่สุด รวมถึง stackoverflow.com/ questions/5963269/ สำหรับคำจำกัดความของ 'ทำซ้ำ' และคำแนะนำเกี่ยวกับวิธีการสร้าง MRE (เป็นไปได้มากว่ามันจะเป็นการง่ายที่สุดในการระบุลิงก์ไปยังชุดข้อมูลทั้งหมดของคุณ ... ) คุณช่วยแสดงด้วยได้ไหม ผลลัพธ์ของ summary(m1) รวมถึงข้อความเตือนใด ๆ ที่อาจถูกพิมพ์ ... ?   -  person Ben Bolker    schedule 12.04.2021
comment
(1) ฉันไม่เห็นว่าคุณจะได้ค่าสัมประสิทธิ์ x3 สำหรับสูง ต่ำ และ reg อย่างไรหากเป็นปัจจัย 3 ระดับ: levels(bd$x3) คืออะไร (2) ในการตรวจสอบข้าม ผลลัพธ์ของ broom.mixed::tidy(bd, effects="fixed", conf.int=TRUE) คืออะไร?   -  person Ben Bolker    schedule 13.04.2021
comment
ฉันรวมบทสรุปของแบบจำลอง ฉันมีข้อผิดพลาด x3 เป็นตัวแปรหมวดหมู่สี่ระดับ ฉันลองใช้โมเดลอื่นที่มีตัวแปรเชิงหมวดหมู่และต่อเนื่องและด้วย ziformula และมันไม่รายงานตัวแปรที่สมบูรณ์ใน zi เช่นกัน อย่างไรก็ตาม ฉันจะพยายามสร้าง MRE สำหรับโมเดลหลายระดับ 3 ระดับ จากนั้นฉันจะอัปเดตโพสต์   -  person MKie45    schedule 13.04.2021
comment
ฉันทำอันหนึ่งสำหรับคุณ   -  person Ben Bolker    schedule 13.04.2021


คำตอบ (1)


tl;dr เท่าที่ฉันสามารถบอกได้ว่านี่เป็นจุดบกพร่องใน confint.glmmTMB (และอาจอยู่ในฟังก์ชันภายใน glmmTMB:::getParms) ในระหว่างนี้ broom.mixed::tidy(m1, effects="fixed") ควรทำสิ่งที่คุณต้องการ (ขณะนี้มีการแก้ไขในเวอร์ชันการพัฒนาบน GitHub ควรทำเป็น CRAN บ้างหรือไม่ เร็ว ๆ นี้ ... )

ตัวอย่างที่ทำซ้ำได้:

ตั้งค่าข้อมูล

set.seed(101)
n <- 1e3
bd <- data.frame(
    year=factor(sample(2002:2018, size=n, replace=TRUE)),
    class=factor(sample(1:20, size=n, replace=TRUE)),
    x1 = rnorm(n),
    x2 = rnorm(n),
    x3 = factor(sample(c("low","reg","high"), size=n, replace=TRUE),
                levels=c("low","reg","high")),
    count = rnbinom(n, mu = 3, size=1))

พอดี

library(glmmTMB)
m1 <- glmmTMB(count~x1+x2+x3+(1|year/class),
          data = bd, zi = ~x2+x3+(1|year/class), family = truncated_nbinom2,
          )

ช่วงความมั่นใจ

confint(m1, "beta_")  ## wrong/ incomplete
broom.mixed::tidy(m1, effects="fixed", conf.int=TRUE)  ## correct

คุณอาจต้องการพิจารณาว่าคุณต้องการช่วงความเชื่อมั่นแบบใด:

  • Wald CI (ค่าเริ่มต้น) คำนวณได้เร็วกว่ามากและโดยทั่วไปก็ใช้ได้ตราบใดที่ (1) ชุดข้อมูลของคุณมีขนาดใหญ่ และ (2) คุณไม่ได้ประมาณค่าพารามิเตอร์ใดๆ ในระดับบันทึก/บันทึกที่อยู่ใกล้กับขอบเขต
  • CIs โปรไฟล์ความน่าจะเป็นมีความแม่นยำมากกว่าแต่ช้ากว่ามาก
person Ben Bolker    schedule 12.04.2021
comment
ฉันลองใช้ broom.mixed แล้ว CI ก็ดูโอเค ขออภัยที่ไม่สามารถทำ MRE ได้ด้วยตัวเองตรงเวลา! แต่ขอบคุณมาก ตอนนี้ฉันเข้าใจมากขึ้นแล้ว ฉันจะอัปโหลดไปยังโพสต์ผลลัพธ์ที่เป็นระเบียบเพื่อตรวจสอบข้าม - person MKie45; 13.04.2021