用户定义的“负指数”链接glm

Question

我试图按照这个例子 在R中修改glm...用户指定的链接函数 进行操作，但是出现了错误。我有二进制数据，并希望将链接函数从“logit”更改为负指数链接。我想预测成功的概率(p) = 1-exp(线性预测器)。

我需要这个链接的原因是在0和0.5之间，p以凸形增加，但“logit”，“cloglog”，“probit”和“cauchy”只允许凹形。请参考附图：预测的p vs 分段观测 模拟数据

location<-as.character(LETTERS[rep(seq(from=1,to=23),30)])
success<-rbinom(n=690, size=1, prob=0.15)
df<-data.frame(location,success)
df$random_var<-rnorm(690,5,3)
df$seedling_size<-abs((0.1+df$success)^(1/df$random_var))
df<-df[order(df$location)]

创建自定义链接函数。注意：eta = 线性预测器，mu = 概率。

negex<-function(){
##link
linkfun<-function(mu) log(-mu+1)
linkinv<-function(eta) 1-exp(eta)
## derivative of inverse link with respect to eta
mu.eta<-function(eta)-exp(eta)
valideta<-function(eta) TRUE
link<-"log(-mu+1)"
structure(list(linkfun=linkfun,linkinv=linkinv,
             mu.eta=mu.eta,valideta=valideta,
             name=link),
        class="link-glm")
}

将种苗大小视为模型成功的函数

negexp<-negex()

model1<-glm(success~seedling_size,family=binomial(link=negexp),data=df)

错误：未找到有效的系数集合，请提供起始值

使用glmer模型（我的最终目标）

model2<-glmer(success~seedling_size+ (1|location),family=binomial(link=negexp),data=df)

错误信息：(function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L, maxit = 100L, maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

我收到不同的错误信息，但我认为无论是使用glmer还是glm，问题都是相同的，那就是我的链接函数出了问题。

- Lake Graboski

Answer 1

我找到了答案。最有帮助的是这个2016年的R讨论串。有两个问题。首先，我的链接函数有误。我将其修改为：

negex <- function() 
 { 
 linkfun <- function(mu) -log(1-mu) 
 linkinv <- function(eta) 1-exp(-eta) 
 mu.eta <- function(eta) exp(-eta) 
 valideta <- function(eta) all(is.finite(eta)&eta>0) 
 link <- paste0("negexp") 
 structure(list(linkfun = linkfun, linkinv = linkinv, 
             mu.eta = mu.eta, valideta = valideta, name = link), 
        class = "link-glm") 
} 

其次，该模型需要特定的起始值。这些值对于您的数据是独一无二的。以下是我找到解决方案时实际使用的数据的前几行：

   site plot sub_plot oak_success oak_o1_gt05ft..1
  0001   10        3           1                0
  0001   12        2           0                0
  0001   12        3           0                0
  0001   12        4           0                0
  0001   13        4           0                0

我不知道如何将完整的数据发布到这个网站，但如果有人想要运行这个示例，请发邮件给我：[e-mail]

（解释：这段文字是一条留言，作者想让读者通过邮件联系他以获得更多信息）

negexp<-negex()

希望这篇文章能对未来的某个人有所帮助，因为我在Stack Overflow或者其他网站上没有找到类似的解决方案。使用新的起始值，我成功地运行了该模型：

starting_values<-c(1,0) #1 for the intercept and 0 for the slope
h_gt05_solo_negex2<-glm(oak_success~ oak_o1_gt05ft..1 , 
                    family=binomial(link=negexp),start=starting_values,data=rocdf)
summary(h_gt05_solo_negex2)
Call:
glm(formula = oak_success ~ oak_o1_gt05ft..1, family = binomial(link = negexp), 
data = lt40, start = starting_values)

Deviance Residuals: 
Min       1Q   Median       3Q      Max  
-1.3808  -0.4174  -0.2637  -0.2637   2.5985  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.034774   0.005484   6.341 2.28e-10 ***
oak_o1_gt05ft..1 0.023253   0.002187  10.635  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1416.9  on 2078  degrees of freedom
Residual deviance: 1213.5  on 2077  degrees of freedom
AIC: 1217.5

Number of Fisher Scoring iterations: 6

收敛存在一些问题。当幼苗高度（oak_o1_gt05ft..1）超过40英尺时，参数估计变得不可靠收敛问题。我在这个范围内观测到非常少的数据，所以我将数据限制在预测值小于40英尺的观测值上，并重新运行模型。我还包括了“站点”（与模拟数据中的“位置”相同）。您在此图中看到的是关于每个站点/位置（黑色圆圈）的橡树成功预测相对于橡树幼苗高度的结果、成功/样本的分组观测值（大绿点）以及没有站点因素的成功概率预测（蓝线）。当考虑站点因素时，看起来幼苗大小变量的斜率更加准确。

不幸的是，我无法在glmer中运行此模型，因此站点必须作为固定效应被纳入，因此，橡树种苗高度的标准误差和斜率估计可能有点保守。

- Lake Graboski