r - 从数字化生存曲线计算危险

Question

我已将已发布的生存曲线数字化并获得以下（生存，时间）点

df<-structure(list(Survival = c(1, 1, 1, 1, 1, 1, 0.99, 0.99, 0.99, 
0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 
0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.98, 0.98, 0.98, 0.98, 
0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 
0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 
0.98, 0.98, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 
0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 
0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 
0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 0.97, 
0.97, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 
0.96, 0.96, 0.96, 0.96, 0.96, 0.96, 0.95, 0.95, 0.95, 0.95, 0.95, 
0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 
0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 
0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 0.94, 
0.94, 0.94, 0.94, 0.94, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 0.93, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 
0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.91, 0.91, 0.91, 0.91, 0.91, 
0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 
0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 
0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 
0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.91, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 
0.9, 0.9, 0.9, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 
0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 0.89, 
0.89, 0.89, 0.89, 0.89, 0.89), time = c(0, 0.34, 0.55, 0.82, 
0.87, 0.87, 0.87, 0.87, 0.87, 1.02, 1.02, 1.24, 1.45, 1.61, 1.74, 
1.83, 1.86, 1.87, 1.87, 1.87, 1.88, 1.96, 2.16, 2.37, 2.75, 3.01, 
3.48, 4.27, 4.5, 4.68, 4.89, 5.19, 5.58, 6.43, 6.75, 7.11, 7.55, 
7.99, 8.41, 8.82, 9.26, 9.72, 10.09, 10.41, 10.68, 10.99, 11.29, 
11.58, 11.85, 12.2, 12.47, 12.69, 13.03, 13.31, 13.72, 14.18, 
14.48, 14.79, 15.26, 15.83, 16.08, 16.55, 16.92, 17.74, 18.54, 
19.24, 19.54, 19.92, 21.14, 21.51, 21.89, 22.98, 23.29, 23.55, 
23.78, 24.1, 24.38, 24.79, 25.09, 25.42, 25.78, 27.09, 27.66, 
28.18, 28.57, 28.97, 29.36, 29.76, 30.15, 30.54, 30.94, 31.33, 
31.73, 32.12, 32.64, 33.19, 33.56, 33.8, 34.05, 34.22, 34.5, 
34.79, 35.1, 35.61, 36.05, 36.47, 36.87, 37.26, 37.65, 38.05, 
38.44, 38.84, 39.23, 39.63, 40.02, 40.42, 40.81, 41.21, 41.6, 
41.99, 42.51, 43.1, 43.41, 43.98, 44.91, 45.55, 46.53, 47.44, 
47.83, 48.17, 49.03, 49.65, 50.45, 50.84, 51.22, 51.61, 52.4, 
52.85, 53.24, 53.62, 54.62, 55.11, 55.46, 56, 56.61, 57.39, 57.75, 
58.15, 58.47, 58.77, 59.01, 59.57, 59.94, 60.44, 60.95, 61.35, 
61.74, 62.14, 62.53, 62.93, 63.32, 63.71, 64.11, 64.62, 65.3, 
65.61, 65.96, 66.21, 66.71, 67.08, 67.6, 68.04, 68.45, 68.85, 
69.24, 69.64, 70.03, 70.43, 70.82, 71.22, 71.61, 72.01, 72.4, 
72.79, 73.19, 73.58, 73.98, 74.37, 74.79, 75.23, 75.68, 76.13, 
76.41, 76.65, 76.88, 77.12, 77.52, 77.88, 78.3, 78.8, 79.22, 
79.64, 79.95, 80.94, 81.32, 81.66, 82.5, 82.86, 83.24, 83.58, 
84.01, 84.27, 84.63, 84.94, 85.28, 86.25, 86.59, 86.91, 88.39, 
88.73, 89.3, 89.59, 90.12, 90.39, 90.67, 91.18, 91.7, 92.14, 
92.55, 92.95, 93.34, 93.74, 94.13, 94.53, 94.92, 95.31, 95.71, 
96.1, 96.5, 96.89, 97.29, 97.68, 98.08, 98.47, 98.87, 99.26, 
99.66, 100.05, 100.44, 100.84, 101.23, 101.63, 102.02, 102.42, 
102.81, 103.21, 103.6, 104, 104.39, 104.78, 105.18, 105.57, 105.97, 
106.36, 106.76, 107.15, 107.55, 107.94, 108.34, 108.73, 109.13, 
109.52, 109.91, 110.31, 110.7, 111.12, 111.56, 112.27, 112.64, 
113.1, 113.59, 114.1, 114.46, 114.8, 115.24, 115.55, 116, 116.34, 
117.03, 117.85, 118.24, 118.58, 118.89, 120.23, 120.55, 120.85, 
121.92, 122.4, 122.89, 123.33, 123.75, 124.14, 124.54, 124.93, 
125.33, 125.72, 126.11, 126.51, 126.9, 127.3, 127.69, 128.09, 
128.48, 128.88, 129.27, 129.67, 130.06, 130.45, 130.85, 131.24, 
131.64, 132.03, 132.43, 132.82, 133.22, 133.61, 134.01, 134.4, 
134.8, 135.19, 135.58, 135.98, 136.37, 136.77, 137.16, 137.56, 
137.95, 138.35, 138.74, 139.14, 139.53, 139.92, 140.34, 140.78, 
141.49, 141.86, 142.29, 142.73, 143.3, 143.81, 144.27, 144.73, 
145.23, 145.66, 146.01, 146.33, 146.92, 147.74, 148.26, 148.66, 
149.21, 149.93, 150.33, 150.87, 151.38, 151.78, 152.17, 152.57, 
152.96, 153.36, 153.75, 154.14, 154.54, 154.93, 155.33, 155.72, 
156.12, 156.51, 156.91, 157.3, 157.81, 158.49, 158.82, 159.33, 
159.8, 160.19, 160.74, 161.25, 161.76, 162.44, 162.76, 163.1, 
163.35, 163.98, 164.38, 164.9, 165.25, 165.57, 166.11, 167.09, 
167.48, 167.9, 168.34, 168.75, 169.15, 169.54, 169.94, 170.35, 
170.79, 171.51, 171.88, 172.31, 172.74, 173.32, 173.82, 174.26, 
174.68, 175.07, 175.47, 175.86, 176.26, 176.65, 177.04, 177.44, 
177.83, 178.23, 178.62, 179.02, 179.41, 179.81, 180.2, 180.6, 
180.99, 181.38, 181.78, 182.17, 182.57, 182.96, 183.36, 183.75, 
184.15, 184.54, 184.94, 185.33, 185.73, 186.12, 186.51, 186.91, 
187.3, 187.7, 188.09, 188.49, 188.88, 189.28, 189.67, 190.07, 
190.46, 190.85, 191.25, 191.64, 192.04, 192.43, 192.83, 193.22, 
193.62, 194.01, 194.41, 194.8, 195.2, 195.59, 195.98, 196.38, 
196.77, 197.17, 197.56, 197.96, 198.37, 198.81, 199.62, 199.91, 
200.33, 200.85, 201.34, 201.84, 202.28, 202.7, 203.09, 203.48, 
203.88, 204.27, 204.67, 205.06, 205.46, 205.85, 206.25, 206.64, 
207.04, 207.43, 207.83, 208.22, 208.61, 209.01, 209.4, 209.8, 
210.19, 210.59, 210.98, 211.38, 211.77, 212.17, 212.56, 212.95, 
213.37, 213.81, 214.33, 214.73, 215.63, 216.28, 216.79, 217.3, 
217.69, 218.09, 218.48, 218.88, 219.27, 219.67, 220.06, 220.46, 
220.85, 221.24, 221.64, 222.03, 222.43, 222.82, 223.22, 223.61, 
224.01, 224.4, 224.8, 225.19, 225.59, 225.98, 226.37, 226.77, 
227.16, 227.56, 227.95, 228.35, 228.74, 229.14, 229.53, 229.93, 
230.32, 230.71, 231.11, 231.5, 231.9, 232.29, 232.81, 233.34, 
233.77, 234.34, 235.23, 235.73, 236.24, 236.64, 237.03, 237.43, 
237.82, 238.22, 238.61, 239.4, 239.77, 240.07, 241.47, 242.05, 
242.56, 242.95, 243.35, 243.74, 244.14, 244.53, 244.93, 245.32, 
245.72, 246.11, 246.51, 246.9, 247.29, 247.69, 248.08, 248.48, 
248.87, 249.27, 249.66, 250.06, 250.45, 250.85, 251.24, 251.63, 
252.03, 252.42, 253.4, 253.99, 254.7, 255.15, 255.56, 256.09, 
256.54, 256.87, 257.19, 258.14, 258.48, 258.8, 259.18, 260.18, 
260.56, 261.05, 261.49, 261.9, 262.3, 262.69, 263.09, 263.48, 
264.14, 264.98, 265.49, 265.87, 266.22, 267.66, 268.15, 268.6, 
269.01, 269.41, 269.8, 270.19, 270.59, 270.98, 271.4, 271.84, 
272.36, 272.76, 273.66, 274.31, 274.82, 275.33, 275.72, 276.12, 
276.63, 277.32, 277.67, 277.91, 278.28, 278.54, 278.84, 279.16, 
279.56, 280, 280.73, 281.48, 282.04, 282.61, 283.31, 284.14, 
284.32, 284.44, 284.64, 284.88, 285.15, 285.45, 285.79, 286.17, 
286.72, 287.15, 287.48, 287.82, 288.35, 289.08, 290.04, 290.32, 
290.82, 291.32, 291.7, 292.2, 292.71, 293.42, 293.94, 294.56, 
295.29, 295.61, 296.27, 297.13, 297.45, 297.79, 298.22, 298.54, 
298.88, 299.41, 299.82, 300.22, 300.61, 301.01, 301.4, 301.8, 
302.19, 302.58, 302.98, 303.37, 303.77, 304.16, 304.56, 305.27, 
306.03, 306.39, 306.86, 307.2, 307.62, 307.95, 308.49, 308.9, 
309.3, 309.69, 310.08, 310.48, 310.87, 311.27, 311.66, 312.06, 
312.45, 312.85, 313.24, 313.64, 314.03, 314.43, 314.82, 315.21, 
315.61, 316, 316.4, 316.79, 317.19, 317.58, 317.98, 318.37, 318.77, 
319.16, 319.55, 319.95, 320.34, 320.74, 321.13, 321.53, 321.92, 
322.32, 323.03, 324.13, 324.49, 325.03, 325.77, 326.1, 326.57, 
326.9, 327.43, 327.85, 328.24, 328.64, 329.03, 329.43, 329.82, 
330.22, 330.61, 331, 331.4, 331.79, 332.19, 332.58, 332.98, 333.37, 
333.79, 334.23, 334.7, 335.03, 335.55, 335.89, 336.36, 336.71, 
337.02, 337.47, 338, 338.51, 338.9, 339.3, 339.69, 340.09, 340.48, 
340.88, 341.27, 341.67, 342.06, 342.45, 342.85, 343.24, 343.64, 
344.03, 344.43, 344.82, 345.22, 345.61, 346.01, 346.4, 346.8, 
347.19, 347.58, 347.98, 348.49, 348.88, 349.24, 349.55, 349.89, 
350.32, 350.65, 350.98, 351.51, 352.25, 353.19, 353.49, 353.83, 
354.27, 354.6, 354.93, 355.46, 355.88, 356.27, 356.67, 357.06, 
357.46, 357.85, 358.25, 358.64, 359.3, 360.16, 360.78, 361.24, 
361.79, 362.27, 362.83, 363.19, 363.79)), class = "data.frame", row.names = c(NA, 
-848L))

我从中得到

plot(df$time, df$Survival, ylim=c(0,1), type="s", lwd=2, bty="n", xlim=c(0,400))

如何从此数据框中计算平滑的瞬时危险函数？我无权访问有风险的号码。

Answer 1

理论上的瞬时危险是累积危险函数的时间导数，它只是 1-Survival。我不同意在这种情况下应该需要“处于风险中的数字”，因为您已经创建了一个经验生存函数。查看生存曲线，您会期望正确的估计有一个早期的峰值，然后逐渐减小。所以看看这个估计器是否足够：

#Pseudo-code
IHF <- { ( S(x+1) - S(x) ) /S(x) }/ # the S(x) in denominator replaces  n-at-risk
            (t(x+1)-t(x)
# R-code
IHF <- function( S, t) { (-diff(S)/head(S, -1))/diff(t) }
plot( head(df$time,-1) , IHF(df$Survival, df$time) )

嗯，并不完全令人满意。你只有 11 个实例有递减，所以我想我明白你为什么想要一个平滑的估计器。让我们看看是否loess足够。没有。无法获得足够的带宽。该muhaz软件包就是为此而设计的。我不确定它是否最适合它，但它是我最熟悉的。它专为审查数据而设计，但其审查状态的默认设置是假设没有案件被审查，因此您只需提供时间，它将考虑生存：

fit2 <- muhaz(head(df$time,-1)[ as.logical(diff(df$Survival)) ], maxtime=360)
plot(fit2)

我对这个结果并不完全满意，因为担心高峰期的到来。所以我尝试了 anotehr 包：

install.packages("bshazard"
library(bshazard
times=head(df$time,-1)[ as.logical(diff(df$Survival)) ]
 times
#[1]   0.87   3.48  13.72  33.80  75.68  90.39 145.66 200.85 272.76 295.29 352.25
  fit<-bshazard( Surv(times,rep(1, length(times))) ~ 1, data=df)

我把那个更好一点。（这有点讽刺，因为我是 pkg:muhaz 的维护者，但我只是一个看守人。）