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

Question

我试着遵循这个示例，修改glm..。用户指定的链接函数，但是我遇到了错误。我有二进制数据，并希望将链接函数从"logit“更改为负指数链接。我想预测

成功概率(P)=1-exp(线性预测器)

我需要这个链接而不是一个内置链接的原因是p在0到0.5之间以凸的方式增加，但是"logit“、"cloglog”、"probit“和"cauchy”只允许凹形。参考见附图：预测p与双值观测

模拟数据

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)

错误在(函数(fr，X，reTrms，家族，nAGQ = 1L，nAGQ= 0L，maxit =100 L，：(maxstephalfit) PIRLS步进半衰期未能减少pwrssUpdate中的偏差)

我得到不同的错误信息，但我认为问题是相同的，无论是使用glmer还是glm，也就是我的链接函数在某种程度上是错误的。

Answer 1

我找到了答案。最有用的是这个R线从2016年开始。有两个问题。首先，我的链接功能是错误的。我将其修改如下：

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

我不知道如何将完整的数据发布到这个站点，但是如果有人想让它运行这个例子，请给我发一封电子邮件: lake.graboski@gmail.com

negexp<-negex()

希望这对将来的人有所帮助，因为我没有发现在堆栈溢出或在线上解决这个问题的其他例子。使用新的起始值，我能够让模型运行：

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英尺的观测值上，并重新运行模型。我还包括“站点”(与模拟数据中的“位置”相同)。您在这个数字中看到的是关于每个地点/位置(黑圈)橡树成功的预测，成功/样本的联合观测(大绿点)和没有场地因素的成功概率的预测(蓝线)。看上去，当考虑到场地因素时，苗木大小变量的斜率更准确。

不幸的是，我无法让这个模型在更好的环境中运行，因此必须将站点作为一个固定的效果，因此，橡树苗高的标准误差和坡度估计可能有点保守。