r - 想要使用 proc nlmixed 将 SAS 代码转换为使用 nlme 的 R 代码

Question

该问题类似于以下帖子中的问题：

“ 麻烦将 proc nlmixed (SAS) 转换为 nlme (R) ”

Answer 1

这看起来根本不像一个混合模型——随机效应是什么？它看起来像一个稍微不寻常的广义线性模型：我认为这可以做到。

glm((1-p)~d1+d2+d3+d4-1,family=binomial(link="log"))

（您必须自己翻转参数的符号）。

蒂瓦里等人。2006 年适合这种类型的模型。

R当量：

dd <- read.table("sas_vals.txt",header=TRUE,na.string=".")
g1 <- glm((1-Y)~d1+d2+d3+d4-1,family=binomial(link="log"),data=dd,
    start=rep(-0.1,4))

这有点困难，可能部分是因为数据集很小（我假设这是真实数据的一个子集——将 4 个参数拟合到 21 个二进制观测值将是一个非常糟糕的主意......）

或者：

library("bbmle")
g2 <- mle2(Y~dbinom(prob=1-exp(-p),size=1),
           parameters=list(p~d1+d2+d3+d4-1),
           start=list(p=0.1),
           data=dd)

对数似然表明mle2拟合更好，glm用这些起始值重新拟合效果更好：

g3 <- update(g1,start=-coef(g2))
g4 <- mle2(Y~dbinom(prob=1-(exp(-p)^exp(b*Treat)),size=1),
           parameters=list(p~d1+d2+d3+d4-1),
           start=list(p=0.1,b=0),
           data=dd)
summary(g4)
##      Estimate Std. Error z value  Pr(z)
## p.d1 0.106163   0.088805  1.1955 0.2319
## p.d2 0.241029   0.178263  1.3521 0.1763
## p.d3 0.105970   0.113455  0.9340 0.3503
## p.d4 0.179232   0.189951  0.9436 0.3454
## b    0.559624   0.740500  0.7557 0.4498

这似乎与 SAS 结果非常吻合。