我正在尝试运行一个
glmer
模型:
这是我正在尝试运行的模型
md$count<-as.integer(md$count)
Model1 <- glmer (count ~ distance_nest_tag + ep_nest + (1|tag_ID),
family = gaussian (link="log"),
data=md )
count
是0到19之间的整数,表示访问次数。distance_nest_tag
是以米为单位的连续变量ep_nest
是二元变量 ( 1, 0)随机效应项只是巢穴和个体的数量。
我也试过这个
Model1 <- glmer (count ~ distance_nest_tag + ep_nest + (1|tag_ID),
family = gaussian (link="log"),
start = coef (lm(md$count ~ md$distance_nest_tag + md$ep_nest)),
data=md)
但我总是明白这个
eval(expr, envir, enclos) 中的错误:找不到有效的起始值:请指定一些
有什么想法吗???
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: EP0 ~ counts + dist_mean_hour + nest_reproductive_period + X.vistors +
(1 | nest_ID) + (1 | tag_ring_ID) + (1 | distance_nesttag_nest) + (1 | date)
Data: only_visits_encounternet
AIC BIC logLik deviance df.resid
34.7643 81.1238 -7.3822 14.7643 752
Random effects:
Groups Name Std.Dev.
date (Intercept) 8.507
distance_nesttag_nest (Intercept) 68.095
tag_ring_ID (Intercept) 64.491
nest_ID (Intercept) 2.251
Number of obs: 762, groups: date, 96; distance_nesttag_nest, 35; tag_ring_ID, 20; nest_ID, 19
Fixed Effects:
(Intercept) counts dist_mean_hour
-13.08952 -1.80786 -0.04528
nest_reproductive_periodI nest_reproductive_periodP X.vistors
-4.04944 -1.56264 0.06601
convergence code 0; 1 optimizer warnings; 0 lme4 warnings
> dput(R)
structure(list(nest_ID = c(6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L), nest_reproductive_period = structure(c(3L,
1L, 2L, 3L, 1L, 2L, 3L, 1L), .Label = c("F", "I", "P"), class = "factor"),
tag_ID = c(23L, 23L, 23L, 4L, 4L, 4L, 5L, 5L), ring_tag_id = structure(c(3L,
3L, 3L, 1L, 1L, 1L, 2L, 2L), .Label = c("AU34180", "AU70442",
"BD12273"), class = "factor"), age_tag = c(5L, 5L, 5L, 4L,
4L, 4L, NA, NA), p3_tag = c(590L, 590L, 590L, 585L, 585L,
585L, 590L, 590L), weight_tag = c(1770L, 1770L, 1770L, 1843L,
1843L, 1843L, 1856L, 1856L), next_box_tag = c(42L, 42L, 42L,
56L, 56L, 56L, 0L, 0L), dist_mean = c(NA, NA, NA, NA, NA,
NA, NA, NA), count = c(3L, 2L, 1L, 0L, 0L, 0L, 0L, 0L), comb = structure(c(2L,
2L, 2L, 3L, 3L, 3L, 1L, 1L), .Label = c("6_0", "6_42", "6_56"
), class = "factor"), distance_nest_tag = c(0.189813484,
0.189813484, 0.189813484, 0.649465717, 0.649465717, 0.649465717,
NA, NA), epp_male_comb = c(NA, NA, NA, NA, NA, NA, NA, NA
), epp_male = c(NA, NA, NA, NA, NA, NA, NA, NA), nest_epmale = c(NA,
NA, NA, NA, NA, NA, NA, NA), neighbour = c(0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L), ep_nest = c(0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L), roost_nextbox_tag = c(8L, 8L, 8L, 335L, 335L, 335L,
332L, 332L), nestroost = structure(c(3L, 3L, 3L, 2L, 2L,
2L, 1L, 1L), .Label = c("6_332", "6_335", "6_8"), class = "factor"),
neighbours_roost = c(1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L)), .Names = c("nest_ID",
"nest_reproductive_period", "tag_ID", "ring_tag_id", "age_tag",
"p3_tag", "weight_tag", "next_box_tag", "dist_mean", "count",
"comb", "distance_nest_tag", "epp_male_comb", "epp_male", "nest_epmale",
"neighbour", "ep_nest", "roost_nextbox_tag", "nestroost", "neighbours_roost"
), class = "data.frame", row.names = c(NA, -8L))
这是一个可重复的小数据集示例,其中响应变量中包含 0 和分组变量。
dat = structure(list(y = c(0, 13.0988072077744, 1.53920266020577, 12.1207857178524,
33.9470080593601, 0, 0, 3.46572339150589, 1.05917038733605, 14.295924854232,
11.2930736701237, 8.1866351980716, 0.0106620794860646, 0.731016675010324,
5.99548577982932, 0, 0, 3.15404516097624, 7.62275500199758, 0.604545763926581,
0, 2.37143378704786, 2.39386320579797, 0.800569675164297), x = c(0,
0, 0, 0, 0, 0, 2.43684896701947, 0.808418724797666, 0.672910050840583,
0, 0.763154948372394, 1.44573423753027, 0.96113385074772, 0,
0.498556550480425, 2.43977373047965, 1.55665618954226, 0.88557694968069,
4.12758995011915, 3.16827587767271, 2.55354765986558, 3.99502024875255,
1.71053826174466, 3.59887218330055), tag = c(1, 1, 1, 1, 1, 1,
2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5)), .Names = c("y",
"x", "tag"), row.names = c(NA, -24L), class = "data.frame")
mod2 = glmer(y ~ x + (1|tag),
data = dat, family = gaussian(link = log),
start = list(theta = 1, coef(lm(y ~ x, data = dat))))
tl;dr 为了使这项工作有效,您需要将
mustart
参数设置为合理的值; mgcv
报告了类似的问题here。在当前版本的 glmer
中,您只需要包含
mustart=pmax(dd$count,1e-3)
在
glmer
通话中。
解释:GLM(迭代重新加权最小二乘法)中通常使用的迭代算法需要设置一个起点,不仅为参数设置起点,还为预测的响应值设置起点。高斯模型的默认设置是将预测值设置为等于观测值;如果观察到的值是特定链接函数的不可能的预测(平均值)值(例如,在对数转换时零值变为负无穷大),那么这将破坏东西(即使观察到的值是合法的(即,即使平均值被限制为正,我们也可以观察到零值)因此,
$initialize
族项目的gaussian(link="log")
分量是:。
if (is.null(etastart) && is.null(start) && is.null(mustart) &&
((family$link == "inverse" && any(y == 0)) || (family$link ==
"log" && any(y <= 0))))
stop("cannot find valid starting values: please specify some")
mustart <- y
所以我们真正想要做的是设置
mustart
,这样所有的值都是正数。
接下来的问题是这是否是正确的统计方法。 这并不疯狂;尝试泊松模型或基于
log(1+x)
的简单线性模型或序数模型也可能是合理的......这取决于您对驱动访问的过程的看法以及您对分布/条件的观察计数分布...