【问题标题】:How to sum over j=1 to (i-1) for each element of [i] (typing formula from article)如何将 [i] 的每个元素的 j=1 与 (i-1) 相加(文章中的键入公式)
【发布时间】:2016-04-14 09:06:20
【问题描述】:

我正在尝试将这个公式输入 R:

公式采用以下输入:

  • M:年死亡人数(全因死亡率);
  • D:每年癌症死亡人数(癌症死亡率);
  • R:每年登记的癌症病例数;
  • N:年中人口规模。
  • w:每个年龄间隔的宽度,例如。 [0-5) 是 5 年宽,最后的间隔是 85+ 年,因此是无限宽的。

上述所有输入向量都是 18 个元素长,因为它们指的是 18 个年龄间隔。 前 17 个年龄区间为 5 年,最后一个区间(85 岁以上)为无限宽。

该公式估计了 Sasieni 等人 2011 年提出的终生癌症风险 http://www.nature.com/bjc/journal/v105/n3/full/bjc2011250a.html

,我不知道怎么打。

下面我尝试在之前和之后实现等式的部分。

# Input data:
M <-   c(140L, 12L, 12L, 59L, 94L, 101L, 117L, 213L, 368L, 607L, 1025L, 
1488L, 2255L, 2787L, 3257L, 3715L, 4231L, 6281L) 


R <-  c(42L, 22L, 28L, 54L, 77L, 108L, 169L, 227L, 293L, 531L, 863L, 
1464L, 2591L, 3334L, 3045L, 2605L, 1890L, 1261L) 


D <-  c(2L, 1L, 2L, 6L, 4L, 7L, 15L, 26L, 67L, 120L, 304L, 497L, 883L, 
1158L, 1321L, 1318L, 1177L, 1065L) 


N <- c(167323L, 168088L, 176017L, 180986L, 168189L, 155506L, 174274L, 
195538L, 207287L, 204711L, 183802L, 174342L, 183415L, 151277L, 
104199L, 71782L, 47503L, 33946L) 

# W width of age interval
w <-  c( 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,Inf )  


# function 
v1 <- numeric()           

for(i in 1:length(R))  {

v1[i] <- R[i] /  ( R[i] + M[i] - D[i] )  *  ( 1 - exp( - (w[i]/N[i]) *  (R[i] + M[i] - D[i]) ) )

}           


sum(v1)

代码看起来尽可能像方程式的答案是首选,以便不了解 R 的同事可以识别代码中的方程式。

答案应该是 0.376127241057822

【问题讨论】:

  • 你能找出问题所在吗?
  • 你知道结果应该是什么?
  • 不,不幸的是我不知道结果应该是什么。
  • 所以当 i=1 时,S_0 是 j=1 到 i=1-1 的总和?或者这被定义为1?或者更可能我不懂数学
  • @rawr 这被定义为零。我刚刚问过math.stackexchange.com/a/1743772/92875

标签: r sum formula


【解决方案1】:

也许这会奏效。论文中没有例子可以查吗?

f <- function(idx) {
  s <- numeric(idx)
  for (i in 1:idx)
    s[i] <- R[i] / (R[i] + M[i] - D[i]) * S(i) * (1 - exp(-w[i] / N[i] * (R[i] + M[i] - D[i])))
  s
}

S <- function(idx) {
  if (idx == 1L)
    return(1)
  s <- numeric(idx - 1)
  for (j in 1:(idx - 1))
    s[j] <- (R[j] + (M[j] - D[j])) / N[j]
  exp(-sum(s))
}

# Input data:
M <-   c(140L, 12L, 12L, 59L, 94L, 101L, 117L, 213L, 368L, 607L, 1025L, 
         1488L, 2255L, 2787L, 3257L, 3715L, 4231L, 6281L) 
R <-  c(42L, 22L, 28L, 54L, 77L, 108L, 169L, 227L, 293L, 531L, 863L, 
        1464L, 2591L, 3334L, 3045L, 2605L, 1890L, 1261L) 
D <-  c(2L, 1L, 2L, 6L, 4L, 7L, 15L, 26L, 67L, 120L, 304L, 497L, 883L, 
        1158L, 1321L, 1318L, 1177L, 1065L) 
N <- c(167323L, 168088L, 176017L, 180986L, 168189L, 155506L, 174274L, 
       195538L, 207287L, 204711L, 183802L, 174342L, 183415L, 151277L, 
       104199L, 71782L, 47503L, 33946L) 
# W width of age interval
w <-  c( 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,Inf )  

f(18)
#  [1] 0.0012516883 0.0006533947 0.0007939380 0.0014874104 0.0022786758 0.0034506651
#  [7] 0.0048088199 0.0057397672 0.0069608906 0.0126706127 0.0226156951 0.0395612334
# [13] 0.0644167605 0.0956951717 0.1184236481 0.1330917708 0.1256574840 0.1421444626

sum(f(18))
# [1] 0.7817021

一个更“R”的方式是

lr <- length(R)
S <- sapply(seq(R), function(idx)
  exp(-sum((R[-(idx:lr)] + (M[-(idx:lr)] - D[-(idx:lr)])) / N[-(idx:lr)])))
sum(R / (R + M - D) * S * (1 - exp(-w / N * (R + M - D))))
# [1] 0.7817021

【讨论】:

  • 这个答案是否考虑到 idx =1 时 s 应该为零? (因此 exp -sum(s))=1 | idx = 1.) 我认为应该......
  • 在 sapply 中使用负索引来获得嵌套总和是一个绝妙的技巧!
【解决方案2】:

也许我读错了问题,但你能通过手动移动 S*0(ai) 向量来解决这个问题吗乘以 1 来说明从 j=1 到 i-1 的总和并与 cumsum 结合?

#df is a data.frame of the example data.  Jump to bottom for code.

#index i = row i
#Using mutate() from dplyr library to make code easier to read
df <- dplyr::mutate(df, RMDN.i = R/(R+M-D) * ( 1 - exp( -(w/N) * (R+M-D) ) ))

#Shift values down one because equation sums from j=1 to i-1.
df$RMDN.i_1 <- c(0, head(df$RMDN.i, -1)) 
df$S0.ai <-exp(-cumsum(df$RMDN.i_1))     #Cumulative sum

#Again, cumulative sum to calculate lifetime risk (Eq. 7)
df <- dplyr::mutate(df, risk = cumsum( R/(R+M-D) * S0.ai * (1 - exp(-(w/N) * (R+M-D)) ) )) 

df
#   age    M    R    D      N   w       RMDN.i     RMDN.i_1     S0.ai        risk
#1    0  140   42    2 167323   5 0.0012516883 0.0000000000 1.0000000 0.001251688
#2    5   12   22    1 168088   5 0.0006540980 0.0012516883 0.9987491 0.001904968
#3   10   12   28    2 176017   5 0.0007949486 0.0006540980 0.9980960 0.002698403
#4   15   59   54    6 180986   5 0.0014896253 0.0007949486 0.9973029 0.004184011
#5   20   94   77    4 168189   5 0.0022834186 0.0014896253 0.9958184 0.006457881
#6   25  101  108    7 155506   5 0.0034612823 0.0022834186 0.9935471 0.009896828
#7   30  117  169   15 174274   5 0.0048298858 0.0034612823 0.9901141 0.014678966
#8   35  213  227   26 195538   5 0.0057738828 0.0048298858 0.9853435 0.020368224
#9   40  368  293   67 207287   5 0.0070171053 0.0057738828 0.9796707 0.027242676
#10  45  607  531  120 204711   5 0.0128095925 0.0070171053 0.9728203 0.039704108
#11  50 1025  863  304 183802   5 0.0229777407 0.0128095925 0.9604383 0.061772810
#12  55 1488 1464  497 174342   5 0.0405424457 0.0229777407 0.9386212 0.099826810
#13  60 2255 2591  883 183415   5 0.0669506082 0.0405424457 0.9013283 0.160171288
#14  65 2787 3334 1158 151277   5 0.1016317397 0.0669506082 0.8429595 0.245842732
#15  70 3257 3045 1321 104199   5 0.1299648254 0.1016317397 0.7614977 0.344810654
#16  75 3715 2605 1318  71782   5 0.1532142188 0.1299648254 0.6686912 0.447263656
#17  80 4231 1890 1177  47503   5 0.1550955224 0.1532142188 0.5737009 0.536242096
#18  85 6281 1261 1065  33946 Inf 0.1946888992 0.1550955224 0.4912792 0.631888708

library(ggplot2)
ggplot(df, aes(x= age, y= risk)) + geom_line() + geom_point() + theme_classic()

# Input data:
df <- data.frame(
        age = seq(0,85, by = 5),  #age band
        M =   c(140L, 12L, 12L, 59L, 94L, 101L, 117L, 213L, 368L, 607L, 1025L, 
                 1488L, 2255L, 2787L, 3257L, 3715L, 4231L, 6281L),
        R =  c(42L, 22L, 28L, 54L, 77L, 108L, 169L, 227L, 293L, 531L, 863L, 
                1464L, 2591L, 3334L, 3045L, 2605L, 1890L, 1261L), 
        D =  c(2L, 1L, 2L, 6L, 4L, 7L, 15L, 26L, 67L, 120L, 304L, 497L, 883L, 
                1158L, 1321L, 1318L, 1177L, 1065L),
        N = c(167323L, 168088L, 176017L, 180986L, 168189L, 155506L, 174274L, 
               195538L, 207287L, 204711L, 183802L, 174342L, 183415L, 151277L, 
               104199L, 71782L, 47503L, 33946L) ,
        w =  c( 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,Inf ) # W width of age interval 
      )

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2019-05-22
    • 2021-07-04
    • 1970-01-01
    • 2021-02-10
    • 1970-01-01
    • 2022-01-13
    • 1970-01-01
    • 2016-03-18
    相关资源
    最近更新 更多