我已经模拟了各种大小和“形状”的对数伽玛数据，然后将伽玛和对数正态模型拟合到这些模拟数据。
这是我的相关代码：
gm_glog <- function(size.i, alpha.i) {
x_i <- runif(size.i, 0, 1) # draw a sample of size 'size'
y.true <- exp(b_0 + b_1*x_i) # produce log gamma data
y_i <- rgamma(size.i, rate = alpha.i/y.true, shape = alpha.i) # random gamma sample
# Gamma Model
log_gamma_model <- glm(y_i ~ x_i, family = Gamma(link = "log"),
control = glm.control(maxit=100, trace = TRUE),
start = c(0.1, 0.2))
log_gamma_summ <- summary(log_gamma_model)
# Lognormal Model
log_norm_model <- glm(y_i ~ x_i, family = gaussian(link = "log"),
control = glm.control(maxit=500, trace = TRUE),
start = c(0.1, 0.2))
log_norm_summ <- summary(log_norm_model)
# DATA FRAME BUILD
data.frame(size = size.i,
alpha = alpha.i,
gamma_mod_int = log_gamma_summ$coefficients["(Intercept)", "Estimate"],
gamma_mod_est = log_gamma_summ$coefficients["x_i", "Estimate"],
gamma_mod_aic = log_gamma_summ$aic,
gamma_mod_dev = log_gamma_summ$deviance.resid[length(log_gamma_summ$deviance.resid)],
gamma_mod_shape = MASS::gamma.shape(log_gamma_model)$alpha,
norm_mod_int = log_norm_summ$coefficients["(Intercept)", "Estimate"],
norm_mod_est = log_norm_summ$coefficients["x_i", "Estimate"],
norm_mod_aic = log_norm_summ$aic,
norm_mod_dev = log_norm_summ$deviance.resid[length(log_norm_summ$deviance.resid)]
)
}
我现在的问题是，我想在一张表中对这些回归结果进行并排比较，其中设计矩阵的每一行[1]对应于函数输出的第一行，再次针对该行[2]，一直到第[40]行。
理想情况下，它看起来像
尺寸 alpha | 摘要gamma glm | 摘要对数正态glm
共有40行，每行分别对应大小和alpha，以最轻松地解释结果。
本质上，我只想合并design.matrix和摘要。
不幸的是，生成一个包含glm摘要的数据帧非常困难，而且我找不到一种像希望的那样逐行合并这些结果的方法。
我已经看到，使用活泼，整洁和一目了然的方式为我提供了所有这些摘要所需的全部信息，但是这两种方式都给我留下了数据帧列表，并且逐行组合它们也使我难以理解。
如果要使用此方法，我仍然想将lapply（model，tidy）的row [1]与lapply（model，glance）的row [1]，lapply（model，tidy）的row [2]与lapply（model，glance）等的row [2]，即使每个列表的行都是不同尺寸的小对象。
我怎样才能最好地做到这一点？有更简单的方法来实现我想要的吗？
编辑：我已经设法获得了单元素列表的偏差残差。仍然不确定如何将这些合并到AIC值等。
考虑使用Map（的包装mapply）的元素明智循环构建数据帧列表，并在每次迭代中运行两个模型，然后将所需的组件提取summary到数据帧：
定义方法
log_models <- function(size.i, alpha.i) {
x_i <- runif(size.i, 0, 1) # draw a sample of size 'size'
y.true <- exp(b_0 + b_1*x_i) # produce log gamma data
y_i <- rgamma(size.i, rate = alpha.i/y.true, shape = alpha.i) # random gamma sample
# Gamma Model
log_gamma_model <- glm(y_i ~ x_i, family = Gamma(link = "log"),
control = glm.control(maxit=100, trace = TRUE),
start = c(0.1, 0.2))
log_gamma_summ <- summary(log_gamma_model)
# Lognormal Model
log_norm_model <- glm(y_i ~ x_i, family = gaussian(link = "log"),
control = glm.control(maxit=500, trace = TRUE),
start = c(0.1, 0.2))
log_norm_summ <- summary(log_norm_model)
# DATA FRAME BUILD
data.frame(size = size.i,
alpha = alpha.i,
gamma_mod_int = log_gamma_summ$coefficients["(Intercept)", "Estimate"],
gamma_mod_est = log_gamma_summ$coefficients["x_i", "Estimate"],
gamma_mod_aic = log_gamma_summ$aic,
gamma_mod_dev = log_gamma_summ$deviance.resid[length(log_gamma_summ$deviance.resid)],
gamma_mod_shape = MASS::gamma.shape(log_gamma_model)$alpha,
norm_mod_int = log_norm_summ$coefficients["(Intercept)", "Estimate"],
norm_mod_est = log_norm_summ$coefficients["x_i", "Estimate"],
norm_mod_aic = log_norm_summ$aic,
norm_mod_dev = log_norm_summ$deviance.resid[length(log_norm_summ$deviance.resid)]
)
}
Map/mapply通话
df_list <- Map(log_models, design.matrix$size, design.matrix$alpha)
# df_list <- mapply(log_models, design.matrix$size, design.matrix$alpha, SIMPLIFY=FALSE)
final_df <- do.call(rbind, df_list)
输出量
final_df
# size alpha gamma_mod_int gamma_mod_est gamma_mod_aic gamma_mod_dev gamma_mod_shape norm_mod_int norm_mod_est norm_mod_aic norm_mod_dev
# 5 5 0.1 -2.39484838 3.808953 2.349387 1.6062347 0.25294152 -0.3943182 0.4366572 21.50163 2.2462398978
# 10 10 0.1 -0.03146698 -1.752435 -48.768787 -2.4685411 0.15839450 -769.8179792 797.7937171 16.72900 0.0073639677
# 15 15 0.1 -6.22434742 11.420125 -146.836144 2.7585789 0.11692945 -0.1601247 1.6135214 102.27202 22.0098432208
# 30 30 0.1 0.26381051 1.067361 -298.873575 -4.7725793 0.08641668 0.2565112 1.0687070 195.59417 -1.7643885736
# 51 5 0.2 -12.23809196 12.760998 -52.109115 0.0412409 0.31666275 -11.1636898 11.2453833 -48.17426 0.0006702163
# 101 10 0.2 1.51817293 -6.261376 -91.417016 -0.7455693 0.12372107 -0.4463434 -1.1394914 31.86825 -0.1580558441
# 151 15 0.2 -0.54878568 3.672312 -17.724359 -1.0910863 0.14922850 -2.7737690 6.2481058 101.48735 0.0621486528
# 301 30 0.2 0.84636917 -1.208503 -25.603596 0.1811917 0.19949756 0.6339933 -0.6533998 168.03056 0.0819567624
# 52 5 0.3 -0.45653740 -2.541001 4.907533 0.8486617 0.66655843 -0.7883221 -0.7289522 10.27774 0.4708082262
# 102 10 0.3 0.70548641 -2.790209 13.450575 0.3375955 0.54226062 1.3245745 -9.0701981 24.19732 -0.8978180162
...
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句