r - 求解 R 中混合模型方程的积分

Question

我正在使用混合效应模型，并且由于我的方法的特殊性我需要解决下面模型的积分，然后制作图表获得的估计值。

换句话说，我需要求解下面的积分:

其中，di^2 是我模型中的 Var3，dh 是混合效应模型对应的函数。

在我所插入的问题的文献中，很少有作品使用用于此目的的混合效应模型，绝大多数仅适用于回归模型简单的线性。但是，对于我的问题，有必要使用混合模型。

模型定义为:

在考虑变量 Var2 的截距中引入了随机效应 bi。

只考虑模型的固定部分，即固定效应模型，我求解积分的过程如下:

数据:https://drive.google.com/file/d/1hFb1OPO0jxQw7_u62swnkRXbOH81ygDD/view?usp=sharing

对于在链接中托管数据，我深表歉意，但是，我找不到包含可能与我的问题匹配的变量的内部 R 数据库。

fitmixedmodel <-  lme(log(Var1)~I(exp(Var3/Var4))+ 
                            (I((Var5/Var4)^3)),
                          random = ~1|Var2, 
                          dados, method="REML")
summary(fitmixedmodel)
volume <- dados[dados$Var5 == 0.1,]
fmixedmodel <- function(Var3, Var5, Var4){
  (pi/40000)*(Var3^2)*(coefficients(summary(fitmixedmodel))[1] + 
                           coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
                           coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3)))
}
vmixedmodel <- function(Var3, Var5, Var4){ 
  integrate(Vectorize(fmixedmodel), lower = 0.1, upper = Var4, Var3 = Var3, Var4 = Var4)$value
  
}
mixed.vol <- mapply(FUN = vmixedmodel,  
                       Var5 = as.list(volume$Var5),
                       Var3 = as.list(volume$Var3),
                       Var4 = as.list(volume$Var4))

所以我得到了下图。

但是请注意，在这个积分的计算中，我没有声明随机效应，也就是说，我只从固定部分积分函数，也没有考虑随机部分。如何解决这个问题，即对调整后的混合模型方程进行实际积分？

Answer 1

我将您的数据集下载为 Data.csv。我必须做一些格式化才能让它在我的本地机器上工作:

library(ggplot2)
library(nlme)
library(data.table)


##################
# Format data   ##
##################


dat <- read.table("Data.csv", 
                  sep=";", 
                  dec=",", 
                  colClasses=c("character",
                                rep("numeric",4)), 
                  skip=1)
setDT(dat)
format(dat,decimal.mark=".")

dat[, Var2 := V1]
dat[, Var3 := as.numeric(V2)]
dat[, Var4 := as.numeric(V3)]
dat[, Var1 := as.numeric(V4)]
dat[, Var5 := as.numeric(V5)]
dat

## this is name used in OP code
dados <- copy(dat[,c("Var2","Var3","Var4","Var1","Var5")])

我稍微重写了代码，以便可以重现您的图形 --这是您的代码，其中有一些小的格式更改:

################### BEGIN OP CODE ####################

fitmixedmodel <-  lme( log(Var1) ~  I(exp(Var3/Var4))+ I((Var5/Var4)^3),
                      random = ~1|Var2, 
                      data = dados, 
                      method="REML")

summary(fitmixedmodel)

volume <- dados[dados$Var5 == 0.1,]

fmixedmodel <- function(Var3, Var5, Var4){
  (pi/40000)*
  (Var3^2)*
  (coefficients(summary(fitmixedmodel))[1] + 
   coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
   coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3)))
}

vmixedmodel <- function(Var3, Var5, Var4){ 
  integrate(Vectorize(fmixedmodel), 
            lower = 0.1, 
            upper = Var4, 
            Var3 = Var3,  
            Var4 = Var4)$value  
}

mixed.vol <- mapply(FUN = vmixedmodel,  
                    Var5 = as.list(volume$Var5),
                    Var3 = as.list(volume$Var3),
                    Var4 = as.list(volume$Var4))

################# END OP CODE ##################

## now verify the graph.  looks good.
ggplot() +
  geom_point(aes(y=mixed.vol, x=volume$Var3, color=volume$Var3))

所以此时我能够重现您的图形。

我最初认为合并随机截距有两种选择。一种是“将它们整合出来”，这将涉及随机截距的方差和二重积分。但事实证明，对于线性回归，这种类型的边缘化不会改变结果。为了向我们自己证明这一点，请看下面的代码，它通过拟合二重积分来积分出服从 Normal(0, 0.1691067^2) 分布的随机截距 b。因为关于 b 的积​​分只能单独隔离 b 并且 E[b] = 0，所以这种方法与 OP 方法没有本质区别。

# Option 1: integrate over the random intercept distribution
# this will require the random intercept variance as well as 
# double integration.

## to be able to accommodate a random intercept, we need to integrate
## over the random intercepts, which are distributed as N(0, sig2)
## where sig2 is 0.1691067^2 as seen from the fitmixedmodel output:

# 
# Random effects:
#   Formula: ~1 | Var2
#          (Intercept)  Residual
# StdDev:   0.1691067 0.2559742

## add "b" random intercept, multiply whole thing by normal density dnorm
integrand <- function(x,  Var3, Var4){
  
  Var5 <- x[1]
  b <- x[2]
  
  (pi/40000)*(Var3^2)*
  (coefficients(summary(fitmixedmodel))[1] + b  + 
   coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
   coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3))) *
    dnorm(b, sd = 0.1691067)
}

vmixedmodel.option1 <- function(Var5, Var3, Var4){ 
  pcubature(integrand, 
            lower = c(0.1,-Inf), 
             upper = c(Var4,Inf), 
             Var3 = Var3, 
             Var4 = Var4)$integral  
}

## this is slow. And unnecessary.  Because the E[b] = 0
mixed.vol.option1 <- mapply(FUN = vmixedmodel.option1,  
                            Var5 = as.list(volume$Var5),
                            Var3 = as.list(volume$Var3),
                            Var4 = as.list(volume$Var4))

max(abs(mixed.vol - mixed.vol.option1))

ggplot() +
  geom_point(aes(y=mixed.vol.option1, x=volume$Var3, color=volume$Var3))

第二种方法是插入估计的随机截距值，很像 OP 方法如何插入 Var4 和 Var3 的值。为了追求这一途径，我们首先创建 volume_ri，它与 volume 数据集相同，但具有 b 的估计值:

## Option 2: plug in the random intercept value.

rand_int <- data.table(Var2 = rownames(fitmixedmodel$coeff$random$Var2), 
                          b = fitmixedmodel$coeff$random$Var2 )
setnames(rand_int, names(rand_int), c("Var2","b"))
rand_int

## merge into `volume` (or `dados` and then re-subset)
volume_ri <- merge(volume,
                   rand_int)

然后基本上我们调整 OP 代码以适应此 b 作为参数或适当的值:

## throw in a b argument
fmixedmodel_ri <- function(Var3, Var5, Var4, b){
  (pi/40000)*
    (Var3^2)*
    (coefficients(summary(fitmixedmodel))[1] + b + 
       coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
       coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3)))
}

## throw in a b argument
vmixedmodel_ri <- function(Var3, Var5, Var4, b){ 
  integrate(Vectorize(fmixedmodel_ri), 
            lower = 0.1, 
            upper = Var4, 
            Var3 = Var3,  
            Var4 = Var4,
               b = b)$value  
}

## plug in the b values
mixed.vol_ri <- mapply(FUN = vmixedmodel_ri,  
                      Var5 = as.list(volume_ri$Var5),
                      Var3 = as.list(volume_ri$Var3),
                      Var4 = as.list(volume_ri$Var4),
                         b = as.list(volume_ri$b))

## now verify the graph. only 8 levels of Var2, so use color
ggplot() +
  geom_point(aes(y=mixed.vol_ri, x=volume_ri$Var3, color=volume_ri$Var2))

---

以下是评论中回答的老问题

老问题:

我担心/困惑的是我评论的行你的意思是Var5 = Var5 - 你介意仔细检查一下吗？并给我留下评论和答案？

此外，还有一个关于随机拦截的单独问题:

是否要对所有随机截距值进行积分

或

是否要为混合模型拟合中的每个唯一 Var2 插入随机截距估计值？