在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))

因此我得到了以下图表。

然而，请注意，计算此积分时我从未声明随机效应，即我仅从固定部分积分函数，也没有考虑随机部分。我该如何解决这个问题，即实际集成调整的混合模型方程？

- user55546

---

有趣的问题，但它不是一个具体的编程问题，因此在这里离题了。我认为stats.stackexchange.com更合适。 - Robert Dodier

嗨，罗伯特，好的！！！ - user55546

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))

所以在这一点上，我能够重现你的图形。

我最初认为有两种选择来纳入随机截距。一种是"将它们整合起来"，这将涉及随机截距的方差和双重积分。但事实证明，在线性回归中，这种边缘化并不会改变结果。为了证明这一点，看看以下代码，它费力地适应一个双重积分来整合随机截距 b，该截距遵循正态分布 Normal(0, 0.1691067^2)。因为对于 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))

---

以下是已回答的旧问题

旧问题：

我的担忧/困惑在于我评论中的那一行，DID YOU MEAN Var5 = Var5 -- 你能再检查一下吗? 并留下一个带有答案的评论吗?

此外，还有一个关于随机截距计算的问题:

你想要在所有随机截距值上进行积分吗？

还是

您想要将混合模型拟合中每个唯一的Var2的随机截距估计插入其中？

- swihart

---

嗨@swihart，感谢您的回复。关于您的问题，即Var3 = Var3是否应该是“Var5 = Var5”，不是的。我会向您解释原因。变量“Var3”指的是控制我正在评估的物体厚度的变量，因此它是“Var3 = Var3”。 参数“upper = Var4”的原因是因为Var4指的是研究对象的总长度。因此，我正在综合考虑每个个体的体积，同时考虑它们的总长度和厚度。因此，我代码中提出的参数是正确的。 - user55546

关于我的问题，就是如何在这个整合中引入随机效应。事实上，我对你描述的两种情况都很感兴趣，即所有随机截距值以及将混合模型拟合的每个唯一Var2的随机截距估计插入其中。所以，只要您介绍这样的解决方案，我就接受答案。感谢您的支持。 - user55546