计算R中的LD50，与SAS中的值不同

ruarlubt 于 12个月前发布在其他

关注(0)|答案(2)|浏览(98)

我在计算昆虫暴露在一种化合物中的LD50，LD90和LD95。我的问题是，我使用不同的程序得到不同的值，我很难找到我的代码是不同的。我使用R中的ecotox包来实现LC_probit函数。
以下是我的数据

Dose Mortality Number
10 12 12
1 8 12
0.5 10 12
0.1 9 12
0.05 11 12
0.01 0 12
0.005 0 12
0.001 0 12

在COVID之前，我可以访问安装了SAS的计算机，这是我使用的代码（包括我的数据）

data CPB;
input conc mort;
number=12;
datalines;
10 12
1 8
0.5 10
0.1 9
0.05 11
0.01 0
0.005 0
0.001 0
0 0
;
*proc sort;
proc probit data=CPB log10 optc inversecl plots=(predpplot ippplot);
model mort/number=conc;
output out=new p=p_hat;
proc print;
run;

下面是我在R中使用的代码（我确保没有在NJ_48中使用控制剂量，因为LC_probit的文档中是这么说的）

NJ_48_probit = LC_probit((Mortality / Number) ~ log10(Dose), data = NJ_48, p = c(50,90,95), weights = Number, het_sig=0.1)

对于LD 50、90和95值，SAS给出0.03092、0.55951和1.27153，而R给出0.0665、0.909和1.91。这种差异从何而来？

来源：https://stackoverflow.com/questions/61873084/calculating-ld50-in-r-different-values-than-in-sas

2条答案

按热度按时间

qij5mzcb1#

在SAS中运行代码会在日志中生成以下警告。
：相对梯度收敛标准0.5340102335大于限值0.0001。这种趋同是值得怀疑的。
注：尽管有上述警告，程序仍在继续。所示结果基于最后一次最大似然迭代。模型拟合的有效性值得怀疑。
第一个警告甚至被复制在结果页面上：Warning from results
这意味着它不能有效地收敛。在那些没有达到收敛的情况下，我认为我们不应该期望SAS和R产生相同的结果。

赞(0）回复(0）举报 12个月前

46scxncf2#

尽管存在收敛误差，但您试图在SAS和R中拟合两个不同的模型。
SAS中的模型包括由optc参数指定的自然响应阈值C。

死亡率/数量= C +（1-C）× Φ（a + b（log 10（剂量）

而R中的不包括自然React阈值C。

死亡率/数量= Φ（a + b（log 10（剂量）

如果不包括自然React阈值C，则R和SAS给予相同的结果。

data CPB;
input conc mort;
number=12;
datalines;
10 12
1 8
0.5 10
0.1 9
0.05 11
0.01 0
0.005 0
0.001 0
0 0
;

proc probit data=CPB log10 inversecl plots=(predpplot ippplot);
model mort/number=conc;
output out=new p=p_hat;
proc print;
run;

| 参数|估计|
| --|--|
| 拦截|一点三二八五|
| Log 10（浓度）|一千二百八十五|

data <- data.frame(conc = c(0, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 10),
           mort = c(0, 0, 0, 0, 11, 9, 10, 8, 12),
           number = 12)
data$conc[data$conc == 0] <- 0.000000001

mod1 <- glm(formula = "mort/number ~ log(conc, 10)", 
            family = binomial(link = "probit"), 
            data = data,
            weights = data$number)

mod1

## Call:  glm(formula = "mort/number ~ log(conc, 10)", family = binomial(link = "probit"), 
##     data = data, weights = data$number)
## 
## Coefficients:
##   (Intercept)  log(conc, 10)  
##         1.329          1.129  
## 
## Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
## Null Deviance:       102.7 
## Residual Deviance: 26.92     AIC: 40.84

coef(mod1)

##   (Intercept) log(conc, 10) 
##      1.328508      1.128505

赞(0）回复(0）举报 12个月前

我来回答

计算R中的LD50，与SAS中的值不同

2条答案

相关问题

热门标签

最新问答