一个联合均值与方差模型的R包——dglm

一、引言

在 R 语言中，dglm 包是用于拟合双参数广义线性模型（Double Generalized Linear Models，简称 DGLMs）的一个工具。这类模型允许同时对均值和方差进行建模，使其适用于处理具有复杂异方差性（方差随均值变化）的数据。dglm 包提供了一个强大的框架，能够扩展传统的广义线性模型（GLM）以包含额外的分布族（双指数族分布，Double Exponential Families and Their Use in Generalized Linear Regression）和链接函数。
特点和应用

• 1.双参数建模：
  dglm 允许用户对响应变量的均值和方差同时进行建模。这一点在数据的条件方差不仅仅是一个常数，而是依赖于其他预测变量时尤为有用。
• 2.支持多种分布：
  包括正态分布、泊松分布、二项分布等，适用于多种类型的数据，如计数数据、比例数据和连续数据等。
  下面我们看一下dglm包如何应用在在联合均值与方差模型中。

二、包的安装与载入

要使用 dglm 包，你首先需要安装并加载它：

install.packages("dglm")
library(dglm)
1
2

下面，我将提供一个实例，演示如何在正态分布假设下使用 dglm 包拟和联合均值和方差模型。这个例子将使用模拟数据，包括响应变量和两个预测变量，来模拟和解释使用过程。
下面的代码生成了响应变量 y，其均值由 x 控制，方差由另一个预测变量 z 控制。

# 安装并加载 dglm 包
if (!requireNamespace("dglm", quietly = TRUE)) {
  install.packages("dglm")
}
library(dglm)
1
2
3
4
5

三、模拟例子

3.1 数据生成

> n <- 100
> p <- 3
> q <- 2
> x <- mvrnorm(n, rep(0, p), diag(rep(1, p)))
> z <- mvrnorm(n, rep(0, q), diag(rep(1, q)))
> Gamma <- rep(1.5, q)
> Beta <- rep(1, p)
> mu <- x %*% Beta
> Sigma <- exp(z %*% Gamma/2)
> y <- rnorm(n, mu, Sigma)
1
2
3
4
5
6
7
8
9
10

3.2 数据查看

> head(x)
           [,1]       [,2]        [,3]
[1,] -1.5780228  0.1444049  0.15079946
[2,] -1.0993938 -0.9157112 -1.21989781
[3,]  0.1035927 -0.3182703  0.51195826
[4,]  0.6475000 -0.1438263 -0.05216163
[5,] -1.0140735 -2.0892163  1.61120512
[6,] -0.4940874 -1.3556227 -2.64037345
> head(z)
           [,1]       [,2]
[1,]  0.1909625 -0.1555610
[2,] -0.2352016  0.9540309
[3,]  0.3996848 -0.3775460
[4,] -0.5754440 -1.5776652
[5,] -1.0378053  1.4773011
[6,]  0.7999732 -0.5783404
> y
  [1]  -1.92185290  -3.71610774   0.23600523   0.27829763   0.79283025
  [6]  -5.83528438   0.75269962  -3.76054139   0.37209131   0.07701730
 [11] -66.38751900   0.08314009   3.28642313 -18.61084616  -2.27184040
 [16]   3.54371296   0.72487268   1.16189009   0.43802927  -4.12017498
 [21]   0.20035594  -0.36651361  -0.75994819  -2.40996426   5.30439360
 [26]  -1.94929237  -1.88685464   1.45781179  -1.17480558  -3.60967067
 [31]  -0.20518913   2.17653752   0.67298245   1.40897920  -0.27651236
 [36]   2.78962790   1.35344365   2.05874102   3.39256784   3.88268947
 [41]  -1.43774535   3.04341943   1.72229299  -0.69458754  -1.59774648
 [46]  -0.11029906  -1.45056626  -2.92491189   0.63004032   1.12543422
 [51]  -0.54342984  -3.19563353   1.23198020   4.19498285  -1.28794485
 [56]   2.74771232   0.30330905  -6.92552297   2.60450584  -2.01859611
 [61]   1.31243200  -0.21181456   0.92944241  -7.21390674  -1.12815181
 [66]  -2.97505852   1.97950406   3.54217820  -0.30837683   1.39522639
 [71]  -3.40358016  -3.52643718  -0.21030846  -5.22556674   1.80455873
 [76]   2.06546172  -2.32681780  -1.51592768   0.39253337  -1.37245595
 [81]  -0.94119131   3.18593841  -5.09087215  -2.48969425   1.25068697
 [86]  -2.26388861  -0.44683619  -3.54009303   1.04009345   1.53777774
 [91]  -0.62336357   0.54738385   1.24649775  -1.28324039   1.35267779
 [96]  -3.84572044  -0.43770251   1.39880786   1.93571598   2.75400275
> 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

3.3 模型估计参数

> library(dglm)
> fit <- dglm(y ~ x[,1] + x[,2] + x[,3] - 1, ~ z[,1] + z[,2] - 1)
> summary(fit)

Call: dglm(formula = y ~ x[, 1] + x[, 2] + x[, 3] - 1, dformula = ~z[, 
    1] + z[, 2] - 1)

Mean Coefficients:
        Estimate Std. Error  t value     Pr(>|t|)
x[, 1] 0.8862072 0.04097214 21.62950 6.881106e-39
x[, 2] 0.8776343 0.03659551 23.98202 1.469803e-42
x[, 3] 1.0764954 0.03109393 34.62076 2.126243e-56
(Dispersion Parameters for gaussian family estimated as below )

    Scaled Null Deviance: 2316.859 on 100 degrees of freedom
Scaled Residual Deviance: 76.88849 on 97 degrees of freedom

Dispersion Coefficients:
       Estimate Std. Error  z value     Pr(>|z|)
z[, 1] 1.761483  0.1413521 12.46167 1.208121e-35
z[, 2] 1.618310  0.1492953 10.83966 2.233068e-27
(Dispersion parameter for Gamma family taken to be 2 )

    Scaled Null Deviance: 522.7625 on 100 degrees of freedom
Scaled Residual Deviance: 125.1396 on 98 degrees of freedom

Minus Twice the Log-Likelihood: 270.2212 
Number of Alternating Iterations: 9 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

其中：输出的 summary(fit) 会提供：

• 均值模型的参数估计、标准误、z 值和 p 值。
• 方差模型的参数估计、标准误、z 值和 p 值。
  使用估计参数进行预测，获取了 β 和 γ 的估计后，你还可以使用这些参数进行预测或进一步的分析。例如，预测新数据点的响应值：
  new_data <- data.frame(x = mvrnorm(n, rep(0, p), diag(rep(1, p)))
  , z = mvrnorm(n, rep(0, q), diag(rep(1, q))))

predicted_values <- predict(fit, newdata = new_data, type = “response”)
print(predicted_values)
这里，predict() 函数使用拟合后的模型和新数据进行预测，type = “response” 表示预测的是响应变量的值。