计算药物LD50用Bliss法最严谨，而改良寇氏法计算的结果误差也不大，因此做了一次改良寇氏法计算LD50的实验。最后需要计算一下结果的可信区间，于是来试试万能的Bootstrap法
安装包（这个例子用不上）
构造样本
组别
剂量 mg/kg
动物数
死亡数
1
110.8
10
0
2
147.7
10
0
3
196.9
10
5
4
262.5
10
8
5
350.0
10
10
某次实验的结果
1
2
3
4
5
6
7
| data <- list(
g1 = rep(0,10),
g2 = rep(0,10),
g3 = c(rep(0,5),rep(1,5)),
g4 = c(rep(0,2),rep(1,8)),
g5 = rep(1,10)
)
|
计算自举置信区间
定义统计量
1
2
3
4
5
6
7
8
9
10
| ln <- log
ld50 <- function(data){
g1 <- mean(sample(x = data$g1, size = 10, replace = T))
g2 <- mean(sample(x = data$g2, size = 10, replace = T))
g3 <- mean(sample(x = data$g3, size = 10, replace = T))
g4 <- mean(sample(x = data$g4, size = 10, replace = T))
g5 <- mean(sample(x = data$g5, size = 10, replace = T))
sigma_p <- sum(g1, g2, g3, g4, g5)
exp(ln(350) - ln(4/3)*(sigma_p - 0.5))
}
|
计算 the bootstrap percentile interval
1
2
3
4
5
6
7
8
9
| set.seed(123)
res <- vector(mode = "list", length = 1000)
for (i in 1:1000){
res[[i]] <- ld50(data)
}
res <- sort(unlist(res))
hist(res)
quantile(res,0.975)
quantile(res,0.025)
|
计算P值
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
| f_Rbisect <- function(lst, value){
low=1
high=length(lst)
if(high == low){return(1)}
mid=length(lst)%/%2
if(lst[low] == value & value == lst[low + 1]){
return(low + 0.5)
}
if(lst[high] == value & value == lst[high - 1]){
return(high - 0.5)
}
if (lst[low] >= value){return(low)}
if (lst[high] <= value){return(high)}
while (lst[mid] != value) {
if (value > lst[mid]){
low <- mid + 1
}else{
high <- mid - 1
}
if (high <= low) { break }
mid <- (low+high)%/%2
}
while(T){
mid0 <- mid - 1
mid2 <- mid + 1
if(lst[mid0] == lst[mid2]){
return(mid)
}
if(lst[mid0] <= value & value <=lst[mid]){
if(lst[mid0] == lst[mid]){
return(mid - 0.5)
}
t = (value - lst[mid0])/(lst[mid] - lst[mid0])
return(mid0 + t)
}
if(lst[mid] <= value & value <= lst[mid2]){
if(lst[mid] == lst[mid2]){
return(mid + 0.5)
}
t = (value - lst[mid])/(lst[mid2] - lst[mid])
return(mid + t)
}
if(value < lst[mid0]){
mid <- mid0
}else{
mid <- mid2
}
}
}
|
1
| f_Rbisect(res, exp(ln(350) - ln(4/3)*(0+0+0.5+0.8+1 - 0.5)))/1000
|