原文链接:http://tecdat.cn/?p=3364
原文出处:拓端数据部落公众号
加载 R 包和数据集
加载包后,咱们将此数据集中蕴含的 12 个情绪变量进行子集化:
mood\_data <- as.matrix(symptom\_data$data\[, 1:12\]) # 子集变量
mood\_labels <- symptom\_data$colnames\[1:12\] # 子集变量标签
colnames(mood\_data) <- mood\_labels
time\_data <- symptom\_data$data_time
对象 mood_data 是一个 1476×12 矩阵,测量了 12 个情绪变量:
> dim(mood_data)
\[1\] 1476 12
> head(mood_data\[,1:7\])
Relaxed Down Irritated Satisfied Lonely Anxious Enthusiastic
\[1,\] 5 -1 1 5 -1 -1 4
\[2,\] 4 0 3 3 0 0 3
\[3,\] 4 0 2 3 0 0 4
\[4,\] 4 0 1 4 0 0 4
\[5,\] 4 0 2 4 0 0 4
\[6,\] 5 0 1 4 0 0 3
time_data 蕴含无关每次测量的工夫戳的信息。数据预处理须要此信息。
> head(time_data)
date dayno beepno beeptime resptime\_s resptime\_e time_norm
1 13/08/12 226 1 08:58 08:58:56 09:00:15 0.000000000
2 14/08/12 227 5 14:32 14:32:09 14:33:25 0.005164874
3 14/08/12 227 6 16:17 16:17:13 16:23:16 0.005470574
4 14/08/12 227 8 18:04 18:04:10 18:06:29 0.005782097
5 14/08/12 227 9 20:57 20:58:23 21:00:18 0.006285774
6 14/08/12 227 10 21:54 21:54:15 21:56:05 0.006451726
该数据集中的一些变量是高度偏斜的,这可能导致不牢靠的参数估计。在这里,咱们通过计算自举置信区间(KS 办法)和可信区间(GAM 办法)来解决这个问题,以判断预计的可靠性。因为本教程的重点是预计时变 VAR 模型,因而咱们不会具体钻研变量的偏度。然而,在实践中,应该在拟合(时变)VAR 模型之前始终查看边际散布。
预计时变 VAR 模型
通过参数 lags = 1,咱们指定拟合滞后 1 VAR 模型,并通过 lambdaSel =“CV”抉择具备穿插验证的参数 λ。最初,应用参数 scale = TRUE,咱们指定在模型拟合之前,所有变量都应标准化。当应用“1 正则化”时,倡议这样做,因为否则参数惩办的强度取决于预测变量的方差。因为穿插验证计划应用随机抽取来定义,因而咱们设置种子以确保重现性。
在查看后果之前,咱们查看了 1476 个工夫点中有多少用于估算,这在调用输入对象的摘要中显示
> tvvar_obj
mgm fit-object
Model class: Time-varying mixed Vector Autoregressive (tv-mVAR) model
Lags: 1
Rows included in VAR design matrix: 876 / 1475 (59.39 %)
Nodes: 12
Estimation points: 20
预计的 VAR 系数的绝对值存储在对象 tvvar_obj $ wadj 中,该对象是维度 p×p×滞后×estpoints 的数组。
参数估计的可靠性
res\_obj <- resample(object = tvvar\_obj,
data = mood_data,
nB = 50,
blocks = 10,seeds = 1:50,
quantiles = c(.05, .95))
res_obj $ bootParameters 蕴含每个参数的教训采样散布。
计算时变预测误差
函数 predict()计算给定 mgm 模型对象的预测和预测误差。
预测存储在 pred\_obj $ 预测中,并且所有时变模型的预测误差组合在 pred\_obj 中:
> pred_obj$errors
Variable Error.RMSE Error.R2
1 Relaxed 0.939 0.155
2 Down 0.825 0.297
3 Irritated 0.942 0.119
4 Satisfied 0.879 0.201
5 Lonely 0.921 0.182
6 Anxious 0.950 0.086
7 Enthusiastic 0.922 0.169
8 Suspicious 0.818 0.247
9 Cheerful 0.889 0.200
10 Guilty 0.928 0.175
11 Doubt 0.871 0.268
12 Strong 0.896 0.195
可视化时变 VAR 模型
可视化下面预计的一部分随工夫变动的 VAR 参数:
# 两个网络图
# 获取均值图的布局
Q <- qgraph(t(mean_wadj), DoNotPlot=TRUE)
saveRDS(Q$layout, "Tutorials/files/layout_mgm.RDS")
#在选定的固定工夫点绘制图形
tpSelect <- c(2, 10, 18)
#
tvvar\_obj$edgecolor\[, , , \]\[tvvar\_obj$edgecolor\[, , , \] == "darkgreen"\] <- c("darkblue")
lty_array <- array(1, dim=c(12, 12, 1, 20))
lty\_array\[tvvar\_obj$edgecolor\[, , , \] != "darkblue"\] <- 2
for(tp in tpSelect) {qgraph(t(tvvar_obj$wadj\[, , 1, tp\]),
layout = Q$layout,
edge.color = t(tvvar_obj$edgecolor\[, , 1, tp\]),
labels = mood_labels,
vsize = 13,
esize = 10,
asize = 10,
mar = rep(5, 4),
minimum = 0,
maximum = .5,
lty = t(lty_array\[, , 1, tp\]),
pie = pred_obj$tverrors\[\[tp\]\]\[, 3\])
}
CIs <- apply(res\_obj$bootParameters\[par\_row\[1\], par_row\[2\], 1, , \], 1, function(x) {quantile(x, probs = c(.05, .95))
} )
# 绘图暗影
polygon(x = c(1:20, 20:1), y = c(CIs\[1,\], rev(CIs\[2,\])), col=alpha(colour = cols\[i\], alpha = .3), border=FALSE)
} #
图 显示了下面预计的时变 VAR 参数的一部分。顶行显示预计点 8,15 和 18 的 VAR 参数的可视化。蓝色实线箭头示意正关系,红色虚线箭头示意负关系。箭头的宽度与相应参数的绝对值成比例。
如果您有任何疑难,请在上面发表评论。
最受欢迎的见解
1.R 语言多元 Logistic 逻辑回归 利用案例
2. 面板平滑转移回归 (PSTR) 剖析案例实现剖析案例实现 ”)
3.matlab 中的偏最小二乘回归(PLSR)和主成分回归(PCR)
4.R 语言泊松 Poisson 回归模型剖析案例
5.R 语言回归中的 Hosmer-Lemeshow 拟合优度测验
6.r 语言中对 LASSO 回归,Ridge 岭回归和 Elastic Net 模型实现
7. 在 R 语言中实现 Logistic 逻辑回归
8.python 用线性回归预测股票价格
9.R 语言如何在生存剖析与 Cox 回归中计算 IDI,NRI 指标