我的任务是评估各种环境变量如何影响年度人口波动。为此，我需要为时间序列计数拟合泊松自回归模型：
其中 N i,j是当年在现场观察到的个体数量，是i当年在现场的环境变量- 这些是输入数据，其余是参数：是一年中现场的预期个体数量，并且是随机的每年的效果。jx_{i,j}ij\mu_{i,j}ij\gamma_{j}
是否可以在 R 中拟合这样的模型？我想避免在贝叶斯框架中拟合它，因为计算需要很长时间（我必须处理 5000 个这样的模型）可能更长。
问问题
5146 次
3 回答
11
首先，让我们创建一些模拟数据（答案末尾的所有临时函数）：
set.seed(12345) # updated to T=20 and L=40 for comparative purposes.
T = 20 # number of years
L = 40 # number of sites
N0 = 100 # average initial pop (to simulate data)
sd_env = 0.8 # to simulate the env (assumed mean 0)
env = matrix(rnorm(T*L, mean=0, sd=sd_env), nrow=T, ncol=L)
# 'real' parameters
alpha = 0.1
beta = 0.05
sd = 0.4
gamma = rnorm(T-1, mean=0, sd=sd)
mu_ini = log(rpois(n=L, lambda=N0)) #initial means
par_real = list(alpha=alpha, beta=beta, gamma=gamma,
sd=sd, mu_ini=mu_ini)
mu = dynamics(par=par_real, x=env, T=T, L=L)
# observed abundances
n = matrix(rpois(length(mu), lambda=mu), nrow=T, ncol=L)
现在，对于一组给定的参数，我们可以模拟每个站点和年份的预期个体数量。由于我们有观察到的个体数量，我们可以为观察写出似然函数（泊松分布），并为增长率的年度偏差添加惩罚（使其呈正态分布）。为此，该函数dynamics将模拟模型并且该函数.getLogLike将计算目标函数。现在我们需要优化目标函数。估计的参数是alpha, beta, 年偏差 ( gamma) 和初始预期个体数 ( mu_ini) , 并且可能sigma.
对于第一次尝试，我们可以为所有参数提供 0 作为初始猜测，但对于我们可以使用初始观察到的丰度的初始预期数字（无论如何这就是 MLE）。
fit0 = fitModel0(obs=n, env=env, T=T, L=L)
Optimal parameters:
alpha beta gamma1 gamma2 gamma3
0.28018842 0.05464360 -0.12904373 -0.15795001 -0.04502903
gamma4 gamma5 gamma6 gamma7 gamma8
0.05045117 0.08435066 0.28864816 0.24111786 -0.80569709
gamma9 gamma10 gamma11 gamma12 gamma13
0.22786951 0.10326086 -0.50096088 -0.08880594 -0.33392310
gamma14 gamma15 gamma16 gamma17 gamma18
0.22664634 -0.47028311 0.11782381 -0.16328820 0.04208037
gamma19 mu_ini1 mu_ini2 mu_ini3 mu_ini4
0.17648808 4.14267523 4.19187205 4.05573114 3.90406443
mu_ini5 mu_ini6 mu_ini7 mu_ini8 mu_ini9
4.08975038 4.17185883 4.03679049 4.23091760 4.04940333
mu_ini10 mu_ini11 mu_ini12 mu_ini13 mu_ini14
4.19355333 4.05543081 4.15598515 4.18266682 4.09095730
mu_ini15 mu_ini16 mu_ini17 mu_ini18 mu_ini19
4.17922360 3.87211968 4.04509178 4.19385641 3.98403521
mu_ini20 mu_ini21 mu_ini22 mu_ini23 mu_ini24
4.08531659 4.19294203 4.29891769 4.21025211 4.16297457
mu_ini25 mu_ini26 mu_ini27 mu_ini28 mu_ini29
4.19265543 4.28925869 4.10752810 4.10957212 4.14953247
mu_ini30 mu_ini31 mu_ini32 mu_ini33 mu_ini34
4.09690570 4.34234547 4.18169575 4.01663339 4.32713905
mu_ini35 mu_ini36 mu_ini37 mu_ini38 mu_ini39
4.08121891 3.98256819 4.08658375 4.05942834 4.06988174
mu_ini40
4.05655031
这可行，但通常一些参数可能是相关的并且更难以从数据中识别，因此可以使用顺序方法（可以阅读Bolker 等人 2013了解更多信息）。在这种情况下，我们逐渐增加参数的数量，从而提高校准每个阶段优化的初始猜测。对于此示例，第一阶段仅估计alpha和beta，并使用对增长率和环境的线性模型的猜测。然后，在第二阶段，我们使用第一次优化的估计值并添加年度偏差作为参数（gamma）。最后，我们使用第二次优化的估计值并将初始期望值作为参数添加。我们最后添加初始期望值，假设初始观察值已经非常接近并且是一个很好的开始猜测，但另一方面，我们不清楚剩余参数的值。
fit = fitModel(obs=n, env=env, T=T, L=L)
Phase 1: alpha and beta only
Optimal parameters:
alpha beta
0.18172961 0.06323379
neg-LogLikelihood: -5023687
Phase 2: alpha, beta and gamma
Optimal parameters:
alpha beta gamma1 gamma2 gamma3
0.20519928 0.06238850 -0.35908716 -0.21453015 -0.05661066
gamma4 gamma5 gamma6 gamma7 gamma8
0.18963170 0.17800563 0.34303170 0.28960181 -0.72374927
gamma9 gamma10 gamma11 gamma12 gamma13
0.28464203 0.16900331 -0.40719047 -0.01292168 -0.25535610
gamma14 gamma15 gamma16 gamma17 gamma18
0.28806711 -0.38924648 0.19224527 -0.07875934 0.10880154
gamma19
0.24518786
neg-LogLikelihood: -5041345
Phase 3: alpha, beta, gamma and mu_ini
Optimal parameters:
alpha beta gamma1 gamma2
0.1962334008 0.0545361273 -0.4298024242 -0.1984379386
gamma3 gamma4 gamma5 gamma6
0.0240318556 0.1909639571 0.1116636126 0.3465693397
gamma7 gamma8 gamma9 gamma10
0.3478695629 -0.7500599493 0.3600551021 0.1361405398
gamma11 gamma12 gamma13 gamma14
-0.3874453347 -0.0005839263 -0.2305008546 0.2819608670
gamma15 gamma16 gamma17 gamma18
-0.3615273177 0.1792020332 -0.0685681922 0.1203006457
gamma19 mu_ini1 mu_ini2 mu_ini3
0.2506129351 4.6639314468 4.7205977429 4.5802529032
mu_ini4 mu_ini5 mu_ini6 mu_ini7
4.4293994068 4.6182382472 4.7039110982 4.5668031666
mu_ini8 mu_ini9 mu_ini10 mu_ini11
4.7610910879 4.5844180026 4.7226353021 4.5823048717
mu_ini12 mu_ini13 mu_ini14 mu_ini15
4.6814189824 4.7130039559 4.6135420745 4.7100006841
mu_ini16 mu_ini17 mu_ini18 mu_ini19
4.4080115751 4.5758092977 4.7209394881 4.5150790425
mu_ini20 mu_ini21 mu_ini22 mu_ini23
4.6171948847 4.7141188899 4.8303375556 4.7392110431
mu_ini24 mu_ini25 mu_ini26 mu_ini27
4.6893526309 4.7237687961 4.8234804183 4.6333012324
mu_ini28 mu_ini29 mu_ini30 mu_ini31
4.6392335265 4.6817044754 4.6260620666 4.8713345071
mu_ini32 mu_ini33 mu_ini34 mu_ini35
4.7107116580 4.5471434622 4.8540773708 4.6129553933
mu_ini36 mu_ini37 mu_ini38 mu_ini39
4.5134108799 4.6231016082 4.5823454113 4.5969785420
mu_ini40
4.5835763300
neg-LogLikelihood: -5047251
比较模型的两个校准，我们可以看到第二个为目标函数提供了较低的值。此外，比较“实际”年度偏差和估计偏差之间的相关性，我们对第二次校准有更高的相关性：
cor(gamma, fit0$par$gamma)
[1] 0.8708379
cor(gamma, fit$par$gamma)
[1] 0.9871758
查看输出，我们可以看到在第一次校准中估计初始期望值（所有站点都被低估）存在一些问题（使用真实数据，通常多阶段校准效果更好）：
par(mfrow=c(3,2), mar=c(3,5,1,1), oma=c(1,1,1,1))
for(i in 1:4) {
ylim=c(0, 1.1*log(max(fit$fitted, n)))
plot(log(fit$fitted[,i]), type="l", col="blue", ylim=ylim,
ylab="mu (log)")
lines(log(fit0$fitted[,i]), col="green")
points(log(mu[,i]), col="red")
mtext(paste("Site ", i), 3, adj=0.05, line=-2)
if(i==3) {
mtext(c("observed", "fitModel0", "fitModel1"), 1, adj=0.95,
line=-1.5:-3.5, col=c("red", "green", "blue"), cex=0.8)
}
}
mus = rbind(mu_ini, fit$par$mu_ini, fit0$par$mu_ini)
barplot(mus, beside=TRUE, col=c("red", "blue", "green"),
ylab="Initial expected population",
xlab="Sites", border=NA)
gam = rbind(gamma, fit$par$gamma, fit0$par$gamma)
barplot(gam, beside=TRUE, col=c("red", "blue", "green"),
ylab="Annual deviates", border=NA)
最后，
system.time(fitModel(obs=n, env=env, T=T, L=L))
user system elapsed
9.85 0.00 9.85
这比@Thierry 使用 INLA（来自 ）提出的解决方案慢了大约四倍summary(model)：
Time used:
Pre-processing Running inla Post-processing Total
0.1070 2.3131 0.0460 2.4661
但是，在对我的函数进行字节编译后，我们得到：
user system elapsed
7.53 0.00 7.53
它快了 24%，现在只比 INLA 方法慢 3 倍。尽管如此，我认为即使对于数千次实验也是合理的（我自己的模型只需要 5 天进行一次优化，所以我在这里可能有偏差）并且由于我们使用的是模拟数据，我们可以比较参数估计的可靠性除了电脑时间。
# The functions -----------------------------------------------------------
require(compiler)
dynamics = function(par, obs, x, T, L) {
alpha = par$alpha
beta = par$beta
gamma = if(!is.null((par$gamma))) par$gamma else rep(0, T-1)
mu_ini = if(!is.null(par$mu_ini)) exp(par$mu_ini) else obs[1,]
mu = matrix(nrow=T, ncol=L)
mu[1,] = mu_ini
for(t in seq_len(T-1)) {
log_mu_new = log(mu[t,]) + alpha + beta*x[t,] + gamma[t]
mu[t+1, ] = exp(log_mu_new)
}
return(mu)
}
dynamics = cmpfun(dynamics)
reListPars = function(par) {
out = list()
out$alpha = as.numeric(par["alpha"])
out$beta = as.numeric(par["beta"])
if(!is.na(par["sd"])) out$sd = as.numeric(par["sd"])
gammas = as.numeric(par[grepl("gamma", names(par))])
if(length(gammas)>0) out$gamma = gammas
mu_inis = as.numeric(par[grepl("mu_ini", names(par))])
if(length(mu_inis)>0) out$mu_ini = mu_inis
return(out)
}
reListPars = cmpfun(reListPars)
.getLogLike = function(par, obs, env, T, L) {
par = reListPars(par)
if(is.null(par$sd)) {
par$sd = if(!is.null(par$gamma)) sd(par$gamma)+0.01 else 1
}
mu = dynamics(par=par, obs=obs, x=env, T=T, L=L)
logLike = sum(obs*log(mu) - mu) - sum(par$gamma^2/(2*par$sd^2))
return(-logLike)
}
.getLogLike = cmpfun(.getLogLike)
.getUpper = function(par) {
par$alpha = 10*par$alpha + 1
par$beta = 10*abs(par$beta) + 1
if(!is.null(par$gamma)) {
if(!is.null(par$sd)) sd = par$sd else sd=sd(par$gamma)
if(sd==0) sd=1
par$gamma = rep(qnorm(0.999, sd=sd), length(par$gamma))
}
if(!is.null(par$mu_ini)) par$mu_ini = 5*par$mu_ini
if(!is.null(par$sd)) par$sd = 10*par$sd
par = unlist(par)
return(par)
}
.getUpper = cmpfun(.getUpper)
.getLower = function(par) {
par$alpha = 0 # alpha>0?
par$beta = -10*abs(par$beta) - 1
if(!is.null(par$gamma)) {
if(!is.null(par$sd)) sd = par$sd else sd=sd(par$gamma)
if(sd==0) sd=1
par$gamma = rep(qnorm(1-0.999, sd=sd), length(par$gamma))
}
if(!is.null(par$mu_ini)) par$mu_ini = 0.2*par$mu_ini
if(!is.null(par$sd)) par$sd = 0.0001*par$sd
par = unlist(par)
return(par)
}
.getLower = cmpfun(.getLower)
fitModel = function(obs, env, T, L) {
r = log(obs[-1,]/obs[-T,])
guess = data.frame(rate=as.numeric(r), env=as.numeric(env[-T,]))
mod1 = lm(rate ~ env, data=guess)
output = list()
output$par = NULL
# Phase 1: alpha an beta only
cat("Phase 1: alpha and beta only\n")
par = list()
par$alpha = as.numeric(coef(mod1)[1])
par$beta = as.numeric(coef(mod1)[2])
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase1 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
# phase 2: alpha, beta and gamma
cat("Phase 2: alpha, beta and gamma\n")
optpar = reListPars(opt$par)
par$alpha = optpar$alpha
par$beta = optpar$beta
par$gamma = rep(0, T-1)
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase2 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
# phase 3: alpha, beta, gamma and mu_ini
cat("Phase 3: alpha, beta, gamma and mu_ini\n")
optpar = reListPars(opt$par)
par$alpha = optpar$alpha
par$beta = optpar$beta
par$gamma = optpar$gamma
par$mu_ini = log(obs[1,])
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par),
control=list(maxit=1000))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase3 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
output$par = reListPars(opt$par)
output$fitted = dynamics(par=output$par, obs=obs, x=env, T=T, L=L)
output$observed = obs
output$env = env
return(output)
}
fitModel = cmpfun(fitModel)
fitModel0 = function(obs, env, T, L) {
output = list()
output$par = NULL
par = list()
par$alpha = 0
par$beta = 0
par$gamma = rep(0, T-1)
par$mu_ini = log(obs[1,])
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase1 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
output$par = reListPars(opt$par)
output$fitted = dynamics(par=output$par, obs=obs, x=env, T=T, L=L)
output$observed = obs
output$env = env
return(output)
}
fitModel0 = cmpfun(fitModel0)
于 2014-03-24T02:49:16.323 回答
1
看看 INLA 包http://www.r-inla.org
它是贝叶斯，但使用集成嵌套拉普拉斯近似，这使得模型的速度可与频率模型（glm，glmm）的速度相媲美。
从Ricardo Oliveros-Ramos 开始mu，envL = 40。首先准备数据集
dataset <- data.frame(
count = rpois(length(mu), lambda = mu),
year = rep(seq_len(T), L),
site = rep(seq_len(L), each = T),
env = as.vector(env)
)
library(reshape2)
n <- as.matrix(dcast(year ~ site, data = dataset, value.var = "count")[, -1])
dataset$year2 <- dataset$year
运行模型
library(INLA)
system.time(
model <- inla(
count ~
env +
f(year, model = "ar1", replicate = site) +
f(year2, model = "iid"),
data = dataset,
family = "poisson"
)
)
user system elapsed
0.18 0.14 3.77
将速度与 Ricardo 的解决方案进行比较
system.time(test <- fitModel(obs=n, env=env, T=T, L=L))
user system elapsed
11.06 0.00 11.06
将速度与常客 glmm 进行比较（无自相关）
library(lme4)
system.time(
m <- glmer(
count ~ env + (1|site) + (1|year),
data = dataset,
family = poisson
)
)
user system elapsed
0.44 0.00 0.44
无自相关的 inla 速度
system.time(
model <- inla(
count ~
env +
f(site, model = "iid") +
f(year, model = "iid"),
data = dataset,
family = "poisson"
)
)
user system elapsed
0.19 0.11 2.09
于 2014-03-24T00:02:36.577 回答
0
模型公式与您给出的不同，但从您的问题标题看来，CRAN 包hhh4中的函数surveillance可能很有趣。它允许您拟合具有随机效应的泊松自回归模型。该函数的文档底部有一些示例。我认为目前的固定效应必须限于每个站点的截距、长期时间趋势和季节性成分，但也许这对你有用。
于 2014-03-22T18:17:56.497 回答