เรียกใช้การถดถอยแบบวนซ้ำสำหรับข้อมูลที่แบ่งออกเป็น N ชิ้นใน R

ฉันมี dataframe ที่มีโครงสร้างดังนี้:

birthwt  tobacco01  pscore  pscoreblocks
3425     0          0.18    (0.177, 0.187]
3527     1          0.15    (0.158, 0.168]
1638     1          0.34    (0.335, 0.345]

คอลัมน์การคลอดบุตรเป็นตัวแปรต่อเนื่องในการวัดน้ำหนักแรกเกิดในหน่วยกรัม คอลัมน์ยาสูบ01มีค่าเป็น 0 หรือ 1 คอลัมน์ pscore มีค่าความน่าจะเป็นระหว่าง 0 ถึง 1 pscoreblocks นำคอลัมน์ pscore และแบ่งออกเป็น 100 บล็อคที่มีขนาดเท่ากัน

ฉันกำลังพยายามหาวิธีที่มีประสิทธิภาพในการทำสิ่งต่อไปนี้สำหรับแต่ละบล็อกใน pscoreblocks ฉันได้รวมโค้ดที่จะใช้งานได้หากฉันเรียกใช้สิ่งนี้กับชุดข้อมูลทั้งหมดโดยไม่ต้องแบ่งพาร์ติชันเป็นบล็อก

1- เรียกใช้การถดถอย

one <- lm(birthwt ~ tobacco01, dfc)

2- หาค่าสัมประสิทธิ์ของตัวแปรยาสูบ01 ในการถดถอย

two <- summary(one)$coefficients[2,1]

3- คูณค่าสัมประสิทธิ์นั้นด้วย: [(จำนวนคนที่ยาสูบ == 1 ในบล็อกนั้น) + (จำนวนคนที่ยาสูบ == 0 ในบล็อกนั้น)] / (จำนวนคนทั้งหมดในกลุ่มนั้น ปิดกั้น)

two_5 <- ((sum(dfc$tobacco01 == 1)) + (sum(dfc$tobacco01 == 0)))/ sum(dfc$tobacco)

three <- two*two_5

4- สุดท้ายนี้ ฉันต้องการที่จะบวกค่าทั้งหมดจาก (3) สำหรับบล็อกทั้งหมด 100 บล็อก

ฉันรู้วิธีทำแต่ละขั้นตอนทีละขั้นตอน แต่ฉันไม่รู้ว่าจะทำซ้ำขั้นตอนเหล่านี้ในบล็อกแยกกันมากกว่า 100 บล็อกได้อย่างไร ฉันลองใช้ group_by(pscoreblocks) แล้วรันการถดถอย แต่ดูเหมือนว่า group_by() และ lm() จะทำงานร่วมกันได้ไม่ดี ฉันได้พิจารณาใช้ pivot_longer() เพื่อสร้างคอลัมน์แยกสำหรับแต่ละบล็อก จากนั้นพยายามเรียกใช้การถดถอยด้วยข้อมูลในรูปแบบนั้น ฉันขอขอบคุณข้อเสนอแนะสำหรับวิธีวนซ้ำทั้งหมด 100 บล็อก

ข้อมูล:

> small <- dput(dfcsmall[1:40,])
structure(list(dbrwt = c(3629, 3005, 3459, 4520, 3095.17811313023, 
3714, 3515, 3232, 3686, 4281, 2645.29691556227, 3714, 3232, 3374, 
3856, 3997, 3515, 3714, 3459, 3232, 3884, 3235, 3008.94507753983, 
3799, 2940, 3389.51332290472, 3090, 1701, 3363, 3033, 2325, 3941, 
3657, 3600, 3005, 4054, 3856, 3402, 2694.09822203382, 3413.03869100037
), tobacco01 = c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 1, 1), pscore = c(0.00988756408875347, 0.183983728674846, 
0.24538311074894, 0.170701594663405, 0.179337494008595,         0.0770304781540708, 
0.164003166666384, 0.0773042518100593, 0.0804603038634144,     0.0611822720382283, 
0.481204657069376, 0.166016137665693, 0.107882394783232,     0.149799473798458, 
0.04130366288307, 0.0360272679038012, 0.476513676221723, 0.214910849480014, 
0.0687582392973688, 0.317662260996216, 0.206183065905609,     0.336553699970873, 
0.0559863953956171, 0.103064791185442, 0.0445362319933672,     0.17097032928289, 
0.245898950803051, 0.146235179401833, 0.284345485401689,     0.152121397241563, 
0.0395696572471225, 0.116669642645446, 0.0672219220193578,     0.297173652687617, 
0.436771917147971, 0.0517299620576624, 0.140760280612358,     0.179726730598874, 
0.0118610298424373, 0.162996197785343), pscoreblocks = structure(c(1L, 
19L, 25L, 18L, 19L, 8L, 17L, 8L, 9L, 7L, 49L, 17L, 11L, 16L, 
5L, 4L, 49L, 22L, 7L, 33L, 21L, 35L, 6L, 11L, 5L, 18L, 25L, 15L, 
29L, 16L, 5L, 12L, 7L, 31L, 45L, 6L, 15L, 19L, 2L, 17L), .Label = c("    [3.88e-05,0.0099]", 
"(0.0099,0.0198]", "(0.0198,0.0296]", "(0.0296,0.0395]", "    (0.0395,0.0493]", 
"(0.0493,0.0592]", "(0.0592,0.069]", "(0.069,0.0789]", "(0.0789,0.0888]", 
"(0.0888,0.0986]", "(0.0986,0.108]", "(0.108,0.118]", "(0.118,0.128]", 
"(0.128,0.138]", "(0.138,0.148]", "(0.148,0.158]", "(0.158,0.168]", 
"(0.168,0.177]", "(0.177,0.187]", "(0.187,0.197]", "(0.197,0.207]", 
"(0.207,0.217]", "(0.217,0.227]", "(0.227,0.237]", "(0.237,0.246]", 
"(0.246,0.256]", "(0.256,0.266]", "(0.266,0.276]", "(0.276,0.286]", 
"(0.286,0.296]", "(0.296,0.306]", "(0.306,0.315]", "(0.315,0.325]", 
"(0.325,0.335]", "(0.335,0.345]", "(0.345,0.355]", "(0.355,0.365]", 
"(0.365,0.375]", "(0.375,0.384]", "(0.384,0.394]", "(0.394,0.404]", 
"(0.404,0.414]", "(0.414,0.424]", "(0.424,0.434]", "(0.434,0.444]", 
"(0.444,0.453]", "(0.453,0.463]", "(0.463,0.473]", "(0.473,0.483]", 
"(0.483,0.493]", "(0.493,0.503]", "(0.503,0.513]", "(0.513,0.522]", 
"(0.522,0.532]", "(0.532,0.542]", "(0.542,0.552]", "(0.552,0.562]", 
"(0.562,0.572]", "(0.572,0.582]", "(0.582,0.591]", "(0.591,0.601]", 
"(0.601,0.611]", "(0.611,0.621]", "(0.621,0.631]", "(0.631,0.641]", 
"(0.641,0.651]", "(0.651,0.66]", "(0.66,0.67]", "(0.67,0.68]", 
"(0.68,0.69]", "(0.69,0.7]", "(0.7,0.71]", "(0.71,0.72]", "(0.72,0.73]", 
"(0.73,0.739]", "(0.739,0.749]", "(0.749,0.759]", "(0.759,0.769]", 
"(0.769,0.779]", "(0.779,0.789]", "(0.789,0.799]", "(0.799,0.808]", 
"(0.808,0.818]", "(0.818,0.828]", "(0.828,0.838]", "(0.838,0.848]", 
"(0.848,0.858]", "(0.858,0.868]", "(0.868,0.877]", "(0.877,0.887]", 
"(0.887,0.897]", "(0.897,0.907]", "(0.907,0.917]", "(0.917,0.927]", 
"(0.927,0.937]", "(0.937,0.946]", "(0.946,0.956]", "(0.956,0.966]", 
"(0.966,0.976]", "(0.976,0.986]"), class = "factor"), blocknumber = c(1L, 
19L, 25L, 18L, 19L, 8L, 17L, 8L, 9L, 7L, 49L, 17L, 11L, 16L, 
5L, 4L, 49L, 22L, 7L, 33L, 21L, 35L, 6L, 11L, 5L, 18L, 25L, 15L, 
29L, 16L, 5L, 12L, 7L, 31L, 45L, 6L, 15L, 19L, 2L, 17L)), row.names =     c(NA, 
-40L), class = c("tbl_df", "tbl", "data.frame"))

r loops apply regression
person melbez    schedule 07.05.2020    source แหล่งที่มา
comment
คุณสามารถรวมรหัสที่คุณใช้สำหรับบล็อก 1 ถึง 4 ทีละรายการได้ไหม   -  person Ronak Shah    schedule 07.05.2020
comment
@RonakShah ใช่ฉันได้เพิ่มมันแล้ว   -  person melbez    schedule 07.05.2020

คำตอบ (2)


arrow_upward
2
arrow_downward

คุณสามารถสร้างฟังก์ชันเพื่อใช้กับ pscoreblocks แต่ละรายการได้

apply_model <- function(data) {
   one <- lm(birthwt ~ tobacco01, data)
   two <- summary(one)$coefficients[2,1]
   two_5 <- ((sum(data$tobacco01 == 1)) + (sum(data$tobacco01 == 0)))/ sum(data$tobacco)
   three <- two*two_5
   return(three)
}

แยกข้อมูลออกเป็น spearate dataframe และใช้ฟังก์ชันนี้กับแต่ละก้อน

library(dplyr)
library(purrr)

dfc %>% group_split(pscoreblocks) %>% map(apply_model)
#OR
#dfc %>% group_split(pscoreblocks) %>% map_dbl(apply_model)

คุณยังสามารถใช้ฐาน R :

lapply(split(dfc, dfc$pscoreblocks), apply_model)

หรือด้วย by :

by(dfc, dfc$pscoreblocks, apply_model)
person Ronak Shah    schedule 07.05.2020
comment
เหตุใดฉันจึงได้รับข้อผิดพลาดต่อไปนี้เมื่อเรียกใช้งานโดยใช้ group_split() และ map() แต่ไม่ได้รับข้อผิดพลาดเมื่อฉันเรียกใช้งานด้วยตัวเอง ข้อผิดพลาดในการสรุป(one)$coefficients[2, 1] : ตัวห้อยอยู่นอกขอบเขต - person melbez; 07.05.2020
comment
ฉันคิดว่าบางรุ่นไม่มี coefficients หรืออย่างน้อย [2,1] มิติไม่ถูกต้อง - person Ronak Shah; 08.05.2020
comment
เมื่อฉันเปรียบเทียบค่าที่ฉันได้รับสำหรับบล็อกแรกเมื่อใช้วิธีนี้ มันแตกต่างจากสิ่งที่ฉันได้รับเมื่อใช้โค้ดบรรทัดเดียวกัน แต่ใช้ตัวกรองเพื่อรวมเฉพาะบล็อกแรกในข้อมูล - person melbez; 12.05.2020
comment
มันจะมีประโยชน์หากคุณเพิ่มข้อมูลบางอย่างโดยใช้ซึ่งฉันสามารถตรวจสอบสิ่งที่คุณพูดได้ - person Ronak Shah; 12.05.2020
comment
ฉันเพิ่มข้อมูลบางส่วน โปรดแจ้งให้เราทราบหากสิ่งนี้ไม่อยู่ในรูปแบบที่เป็นประโยชน์ เนื่องจากฉันไม่เคยทำเช่นนี้มาก่อน ฉันใช้ dput() เพื่อรับสิ่งนี้ - person melbez; 12.05.2020
comment
มีคนตอบกลับที่นี่และทราบว่าเกิดอะไรขึ้น: stackoverflow.com/questions/61741207/ - person melbez; 12.05.2020
comment
@melbez ดังนั้นฟังก์ชั่นทำงานได้ตามที่คาดไว้ใช่ไหม? มันเป็นเพียงเรื่องของการตีความ - person Ronak Shah; 12.05.2020

arrow_upward
2
arrow_downward

คำถามน่าจะเป็นโมดูลโครงการ

ฉันเชื่อว่าปัญหาหลักสองประการในคำถามคือ 1 และ 2 ดังนั้นการตอบคำถามเหล่านั้น

ขั้นตอน:

  1. ซ้อนชุดข้อมูลของคุณโดยใช้ pscoreblocks

    d_nested <- d %>% group_by(pscoreblocks) %>% nest()

  2. เขียนฟังก์ชันให้กับโมเดล

    mod_fun <- function(df){ lm( birthwt ~ tobacco01, data = df) }

  3. ใช้ฟังก์ชันด้านบนเพื่อสร้างแบบจำลอง

    m_d <- d_nested %>% mutate(model = map(data, mod_fun))

  4. สร้างฟังก์ชันอื่นเพื่อแยกค่าสัมประสิทธิ์ของแต่ละรุ่น

    b_fun <- function(mod){ coefficients(mod)[[1]] }

  5. สุดท้ายให้ใช้ฟังก์ชันข้างต้น

    m_d %>% transmute(coeff = map_dbl(model, b_fun))

จะให้ผลลัพธ์แก่คุณ [coeffs เหมือนกับ data เพราะเรามีจุดข้อมูลเพียงจุดเดียวต่อกลุ่ม] เป็น

# A tibble: 3 x 2
# Groups:   pscoreblocks [3]
  pscoreblocks   coeff
  <chr>          <dbl>
1 (0.177, 0.187]  3425
2 (0.158, 0.168]  3527
3 (0.335, 0.345]  1638

ข้อมูล:

structure(list(birthwt = c(3425, 3527, 1638), tobacco01 = c(0, 
1, 1), pscore = c(0.18, 0.15, 0.34), pscoreblocks = c("(0.177, 0.187]", 
"(0.158, 0.168]", "(0.335, 0.345]")), row.names = c(NA, -3L), class = c("tbl_df", 
"tbl", "data.frame")) -> d
person nikn8    schedule 07.05.2020