目录
- 一、引言
- 二、包的安装与载入
- 三、模拟例子
- 3.1 数据生成
- 3.2 数据查看
- 3.3 模型估计参数
一、引言
在 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)
下面，我将提供一个实例，演示如何在正态分布假设下使用 dglm 包拟和联合均值和方差模型。这个例子将使用模拟数据，包括响应变量和两个预测变量，来模拟和解释使用过程。
下面的代码生成了响应变量 y，其均值由 x 控制，方差由另一个预测变量 z 控制。
# 安装并加载 dglm 包
if (!requireNamespace("dglm", quietly = TRUE)) {
install.packages("dglm")
}
library(dglm)
三、模拟例子
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)
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
>
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
其中：输出的 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” 表示预测的是响应变量的值。
