解决R中混合模型方程的积分

Question

我正在使用混合效果模型，并且由于我的方法的细节 我需要求解下面模型的积分，然后进行图形 获得的估计值。

换句话说，我需要解决以下积分：

其中， di^2 是 var3 在我的模型中， dh 是与混合效果相对应的功能模型。

在有关我插入的问题的文献中，几乎没有作品 为此目的的混合效应模型，绝大多数仅与回归模型一起工作 简单线性。但是，就我的问题而言，有必要使用混合模型。

该模型由：

考虑变量 var2 ，在截距中引入了随机效果。

仅考虑模型的固定部分，即是固定效应模型，我执行的解决积分的过程如下：

data： https：//drive.google.com/file.com/file.com/file/d/1hfb1hfb1hfb1hfb1opo链接

，但是，我找不到具有可能与我的问题相匹配的变量的内部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))

因此，此时我能够复制您的图形。

我最初认为有两个选择随机拦截。一种是“整合它们”，这将涉及随机截距的方差和双重积分。但是事实证明，对于线性回归，这种边缘化不会改变结果。为了向我们证明这一点，请查看以下代码，这些代码遇到了拟合双重积分以集成遵循正常（0，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))

第二种方法是插入估计的随机拦截值，就像 var4 和 var3 中的OP方法插入值一样。要追求这一途径，我们首先创建 volume_ri ，它与卷数据集相同，但具有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 = var5 - 您是否介意仔细检查一下。并给我回答的评论吗？

另外，关于考虑随机截距的另一个问题：

您是否要在所有随机拦截值

或

上集成，您是否想插入混合模型拟合中每个唯一的VAR2的随机截距估计？