r - 使用 gamma 系列测试 glm 时出错

Question

我目前正在研究 10 年内六个法国海湾的贝类毒性。我创建了一个代表这一年出现的毒素数量的代理。
现在我想通过各种参数来解释这种毒性。目前，我专注于研究微藻大量繁殖对毒性的影响。为了做到这一点，我描述了每年一次的盛开:一年中观察到的最大细胞数、盛开的持续时间、开始的一周等......

因为我的毒性数据很少，所以我想让事情变得非常简单，所以我在这些绽放参数中的 1、2 或 3 个参数上测试 glms，并且我允许交互。

该代码适用于对数正态族，但我想与 Gamma 分布和 进行比较。甚至是准 Gamma (均值方差在 var=2*mean² 的图形上效果很好)。

当我从 gamma 开始时，我收到以下消息:

 Error: no valid set of coefficients has been found: please supply starting values
In addition: Warning message:
In log(ifelse(y == 0, 1, y/mu)) : NaNs produced

它出现在第七个模型中:一个带有 的模型只有一个 参数:代理~DureeB!

我读到我需要用 start= 找到 coef参数，但我有 15 000 个模型要检查，包含 1、2 或 3个参数。所以这里是我的问题:

为什么在这么简单的模型上不起作用？？

是否有我遗漏的元素，我的代码中的错误解释了该消息？

有没有办法为我的所有模型设置起始系数？

这是我的代码和数据:

output<-sapply(models,  function (x) { # x<-models[7]

    glm.01<-glm(as.formula(x), family=Gamma, data=TestGLM)

    Aic<-AIC(glm.01)  #+2*sum(log(TestGLM$Proxy))# correction for log normal AIC
    Dev<-(glm.01$null.deviance-glm.01$deviance)/(glm.01$null.deviance)
    c(Aic, Dev)

  })

> dput(TestGLM)

structure(list(Annee = c(2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 
2009L, 2010L, 2011L, 2012L, 2003L, 2004L, 2005L, 2006L, 2007L, 
2008L, 2009L, 2010L, 2011L, 2012L, 2003L, 2004L, 2005L, 2006L, 
2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2003L, 2004L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2003L, 2004L, 
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2003L, 
2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L
), Proxy = c(1, 1, 4.7, 72.6, 37.3, 32.5, 12.4, 164.2, 33.3, 
63.3, 0.85, 10.6, 37.7, 21.9, 37.4, 30.2, 38.3, 445.9, 57, 29.7, 
0.48, 48.4, 31.7, 0.85, 2.2, 0.85, 2.1, 1.9, 75.1, 287.9, 0.15, 
0.85, 0.85, 2.65, 1, 0.85, 1, 0.05, 1, 0.45, 0.15, 1, 1, 0.85, 
0.85, 0.85, 0.1, 194.85, 26.1, 21.8, 1, 46.5, 12, 4, 180.6, 31.7, 
11.7, 3.2, 6, 0.7), Libellé = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("BaiedeConcarneau", "BaiedeQuiberon", 
"BaiedeSeine", "BaiedeStBrieuc", "PertuisBreton", "RadedeBrest"
), class = "factor"), Max = c(462300, 1616150, 1378900, 9362300, 
718900, 5521866.667, 742400, 258550, 486700, 1360400, 58700, 
1200000, 2e+06, 255300, 252000, 614500, 1263900, 383100, 440833.3333, 
736700, 1173300, 218950, 971300, 468566.6667, 176300, 488400, 
696000, 28850, 65700, 575650, 11000, 8900, 13100, 36200, 133500, 
111300, 5500, 5500, 8900, 5800, 33800, 51500, 19100, 96700, 13800, 
13000, 4900, 450766.6667, 4200, 30700, 47900, 393950, 33000, 
62700, 325700, 1304000, 290100, 116600, 85600, 62300), QuinzaineMax = c(13L, 
13L, 7L, 9L, 19L, 13L, 12L, 9L, 15L, 13L, 12L, 12L, 7L, 12L, 
12L, 12L, 12L, 9L, 11L, 13L, 12L, 11L, 12L, 16L, 14L, 11L, 10L, 
15L, 17L, 11L, 12L, 16L, 11L, 20L, 11L, 11L, 14L, 15L, 9L, 6L, 
14L, 12L, 14L, 10L, 21L, 12L, 8L, 7L, 9L, 9L, 12L, 12L, 15L, 
15L, 19L, 10L, 12L, 17L, 10L, 15L), DébutB = c(11L, 13L, 7L, 
9L, 11L, 13L, 10L, 8L, 15L, 7L, 11L, 12L, 7L, 10L, 11L, 10L, 
12L, 9L, 11L, 7L, 12L, 11L, 11L, 15L, 13L, 10L, 9L, 4L, 7L, 10L, 
9L, 7L, 7L, 4L, 10L, 10L, 4L, 3L, 1L, 6L, 9L, 11L, 12L, 6L, 6L, 
1L, 1L, 6L, 1L, 8L, 11L, 12L, 7L, 10L, 9L, 10L, 8L, 10L, 10L, 
8L), FinB = c(14L, 15L, 15L, 9L, 19L, 13L, 18L, 26L, 15L, 14L, 
17L, 13L, 13L, 17L, 18L, 13L, 12L, 17L, 12L, 13L, 12L, 13L, 13L, 
16L, 21L, 11L, 10L, 16L, 18L, 17L, 22L, 17L, 20L, 20L, 11L, 11L, 
17L, 26L, 22L, 26L, 15L, 13L, 17L, 14L, 26L, 18L, 21L, 10L, 22L, 
12L, 16L, 19L, 21L, 18L, 19L, 10L, 14L, 20L, 18L, 18L), DuréeB = c(1L, 
1L, 1L, 1L, 4L, 1L, 3L, 3L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 4L, 1L, 
2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 2L, 1L, 1L, 3L, 2L, 2L, 1L, 
1L, 2L, 2L, 5L, 2L, 3L, 4L, 1L, 3L, 4L, 1L, 1L, 1L, 3L, 2L, 1L, 
5L, 6L, 1L, 1L, 1L, 1L, 1L, 1L, 8L, 1L, 1L), AbCum = c(1232834.481, 
2295533.333, 2618400, 9540123.181, 2721369.612, 6402845.628, 
2131316.116, 1098695.022, 635295.5066, 2341036.687, 176743.6525, 
2008022.948, 4088103.498, 691735.877, 564020.5998, 1396076.209, 
1701772.35, 699439.2121, 667583.3333, 1274466.667, 1222289.094, 
501194.8149, 1882821.151, 645556.1372, 317050, 641859.2929, 1094461.592, 
102315.0749, 252271.4604, 1719625.936, 27751.20166, 23800, 56531.29739, 
90193.6705, 189335.9309, 178272.214, 15000, 26350, 47748.27689, 
25566.68521, 78063.53164, 97687.20404, 69380.95845, 188604.3622, 
48992.92008, 32453.63676, 29969.77734, 1003741.845, 20519.75684, 
83103.25321, 168584.9387, 481250.7282, 92726.1788, 128114.5066, 
811650, 1526211.659, 546353.7302, 440655.1676, 232500, 219290.8449
), NbBloom = c(2L, 2L, 3L, 1L, 2L, 1L, 3L, 4L, 1L, 3L, 3L, 1L, 
3L, 3L, 2L, 1L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 3L, 
4L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 2L, 4L, 5L, 5L, 2L, 1L, 2L, 3L, 
6L, 3L, 3L, 2L, 7L, 1L, 1L, 2L, 3L, 3L, 3L, 1L, 2L, 2L, 4L, 2L
)), .Names = c("Annee", "Proxy", "Libellé", "Max", "QuinzaineMax", 
"DébutB", "FinB", "DuréeB", "AbCum", "NbBloom"), class = "data.frame", row.names = c("4", 
"10", "16", "22", "28", "34", "40", "46", "52", "58", "5", "11", 
"17", "23", "29", "35", "41", "47", "53", "59", "1", "7", "13", 
"19", "25", "31", "37", "43", "49", "55", "2", "8", "14", "20", 
"26", "32", "38", "44", "50", "56", "6", "12", "18", "24", "30", 
"36", "42", "48", "54", "60", "3", "9", "15", "21", "27", "33", 
"39", "45", "51", "57"))

我想测试的一些模型

models<-c("Proxy~ Libellé", "Proxy~ Annee", "Proxy~ Max", "Proxy~ QuinzaineMax", 
"Proxy~ DébutB", "Proxy~ FinB", "Proxy~ DuréeB", "Proxy~ AbCum", 
"Proxy~ NbBloom", "Proxy~ Libellé:Annee", "Proxy~ Libellé:Max", 
"Proxy~ Libellé:QuinzaineMax", "Proxy~ Libellé:DébutB", "Proxy~ Libellé:FinB", 
"Proxy~ Libellé:DuréeB", "Proxy~ Libellé:AbCum", "Proxy~ Libellé:NbBloom", 
"Proxy~ Annee:Max", "Proxy~ Annee:QuinzaineMax", "Proxy~ Annee:DébutB", 
"Proxy~ Annee:FinB", "Proxy~ Annee:DuréeB", "Proxy~ Annee:AbCum", 
"Proxy~ Annee:NbBloom", "Proxy~ Max:QuinzaineMax", "Proxy~ Max:DébutB", 
"Proxy~ Max:FinB", "Proxy~ Max:DuréeB", "Proxy~ Max:AbCum", "Proxy~ Max:NbBloom", 
"Proxy~ QuinzaineMax:DébutB", "Proxy~ QuinzaineMax:FinB", "Proxy~ QuinzaineMax:DuréeB", 
"Proxy~ Libellé+Libellé:Max", "Proxy~ Libellé+Libellé:QuinzaineMax", 
"Proxy~ Libellé+Libellé:DébutB", "Proxy~ Libellé+Libellé:FinB", 
"Proxy~ Libellé+Libellé:DuréeB", "Proxy~ Libellé+Libellé:AbCum", 
"Proxy~ Libellé+Libellé:NbBloom", "Proxy~ Libellé+Annee:Max", 
"Proxy~ Libellé+Annee:QuinzaineMax", "Proxy~ Libellé+Annee:DébutB", 
"Proxy~ Libellé+Annee:FinB", "Proxy~ Libellé+Annee:DuréeB", "Proxy~ Libellé+Annee:AbCum", 
"Proxy~ Libellé+Annee:NbBloom", "Proxy~ Libellé+Max:QuinzaineMax", 
"Proxy~ Libellé+Max:DébutB", "Proxy~ Libellé+Max:FinB", "Proxy~ Libellé+Max:DuréeB", 
"Proxy~ Libellé+Max:AbCum", "Proxy~ Libellé+Max:NbBloom", "Proxy~ Libellé+QuinzaineMax:DébutB", 
"Proxy~ Libellé+QuinzaineMax:FinB", "Proxy~ Libellé+QuinzaineMax:DuréeB", 
"Proxy~ Libellé+QuinzaineMax:AbCum", "Proxy~ Libellé+QuinzaineMax:NbBloom", 
"Proxy~ Libellé+DébutB:FinB", "Proxy~ Libellé+DébutB:DuréeB", 
"Proxy~ Libellé+DébutB:AbCum", "Proxy~ Libellé+DébutB:NbBloom", 
"Proxy~ Libellé+FinB:DuréeB", "Proxy~ Libellé+FinB:AbCum", "Proxy~ Libellé+FinB:NbBloom", 
"Proxy~ Libellé+DuréeB:AbCum", "Proxy~ Libellé+DuréeB:NbBloom", 
"Proxy~ Libellé+AbCum:NbBloom", "Proxy~ Annee+Max", "Proxy~ Annee+QuinzaineMax", 
"Proxy~ Annee+DébutB", "Proxy~ Annee+FinB", "Proxy~ Annee+DuréeB", 
"Proxy~ Annee+AbCum", "Proxy~ Annee+NbBloom", "Proxy~ Annee+Libellé:Annee", 
"Proxy~ Annee+Libellé:Max", "Proxy~ Annee+Libellé:QuinzaineMax", 
"Proxy~ Annee+Libellé:DébutB", "Proxy~ Annee+Libellé:FinB", "Proxy~ Annee+Libellé:DuréeB", 
"Proxy~ Annee+Libellé:AbCum", "Proxy~ Annee+Libellé:NbBloom", 
"Proxy~ Annee+Annee:Max", "Proxy~ Annee+Annee:QuinzaineMax", 
"Proxy~ Annee+Annee:DébutB", "Proxy~ Annee+Annee:FinB", "Proxy~ Annee+Annee:DuréeB", 
"Proxy~ Annee+Annee:AbCum", "Proxy~ Annee+Annee:NbBloom", "Proxy~ Annee+Max:QuinzaineMax", 
"Proxy~ Annee+Max:DébutB", "Proxy~ Annee+Max:FinB", "Proxy~ Annee+Max:DuréeB", 
"Proxy~ Annee+Max:AbCum", "Proxy~ Annee+Max:NbBloom", "Proxy~ Annee+QuinzaineMax:DébutB", 
"Proxy~ Libellé+Libellé:DuréeB+Max:DébutB", "Proxy~ Libellé+Libellé:DuréeB+Max:FinB", 
"Proxy~ Libellé+Libellé:DuréeB+Max:DuréeB", "Proxy~ Libellé+Libellé:DuréeB+Max:AbCum", 
"Proxy~ Libellé+Libellé:DuréeB+Max:NbBloom", "Proxy~ Libellé+Libellé:DuréeB+QuinzaineMax:DébutB", 
"Proxy~ Libellé+Libellé:DuréeB+QuinzaineMax:FinB", "Proxy~ Libellé+Libellé:DuréeB+QuinzaineMax:DuréeB", 
"Proxy~ Libellé+Libellé:DuréeB+QuinzaineMax:AbCum", "Proxy~ Libellé+Libellé:DuréeB+QuinzaineMax:NbBloom", 
"Proxy~ Libellé+Libellé:DuréeB+DébutB:FinB", "Proxy~ Libellé+Libellé:DuréeB+DébutB:DuréeB", 
"Proxy~ Libellé+Libellé:DuréeB+DébutB:AbCum", "Proxy~ Libellé+Libellé:DuréeB+DébutB:NbBloom", 
"Proxy~ Libellé+Libellé:DuréeB+FinB:DuréeB", "Proxy~ Libellé+Libellé:DuréeB+FinB:AbCum", 
"Proxy~ Libellé+Libellé:DuréeB+FinB:NbBloom", "Proxy~ Libellé+Libellé:DuréeB+DuréeB:AbCum", 
"Proxy~ Libellé+Libellé:DuréeB+DuréeB:NbBloom", "Proxy~ Libellé+Libellé:DuréeB+AbCum:NbBloom", 
"Proxy~ Libellé+Libellé:AbCum+Libellé:NbBloom", "Proxy~ Libellé+Libellé:AbCum+Annee:Max", 
"Proxy~ Libellé+Libellé:AbCum+Annee:QuinzaineMax", "Proxy~ Libellé+Libellé:AbCum+Annee:DébutB", 
"Proxy~ Libellé+Max:QuinzaineMax+QuinzaineMax:DébutB", "Proxy~ Libellé+Max:QuinzaineMax+QuinzaineMax:FinB", 
"Proxy~ Libellé+Max:QuinzaineMax+QuinzaineMax:DuréeB", "Proxy~ Libellé+Max:QuinzaineMax+QuinzaineMax:AbCum", 
"Proxy~ Libellé+Max:QuinzaineMax+QuinzaineMax:NbBloom", "Proxy~ Libellé+Max:QuinzaineMax+DébutB:FinB", 
"Proxy~ Libellé+Max:QuinzaineMax+DébutB:DuréeB", "Proxy~ Libellé+Max:QuinzaineMax+DébutB:AbCum", 
"Proxy~ Libellé+Max:QuinzaineMax+DébutB:NbBloom", "Proxy~ Libellé+Max:QuinzaineMax+FinB:DuréeB", 
"Proxy~ Libellé+Max:QuinzaineMax+FinB:AbCum", "Proxy~ Libellé+Max:QuinzaineMax+FinB:NbBloom", 
"Proxy~ Libellé+Max:QuinzaineMax+DuréeB:AbCum", "Proxy~ Libellé+Max:QuinzaineMax+DuréeB:NbBloom", 
"Proxy~ Libellé+Max:QuinzaineMax+AbCum:NbBloom", "Proxy~ Libellé+Max:DébutB+Max:FinB") 

感谢您的时间!

Answer 1

注意 Gamma 的默认链接函数家庭是反向链接。如果您要与对数正态模型进行比较，如果我理解正确，这意味着对响应进行对数转换并使用线性回归，您可能需要一个对数链接( family=Gamma(link=log) )。

根据我的经验，日志链接也往往不那么脆弱，因此这可以很好地解决您的收敛问题。事实上，当我使用对数链接拟合模型时，它按预期收敛。

它不会解决您试图使 15000 个模型变得合理的其他问题，但这是另一个问题。