I'm running a generalised additive mixed model using the gamlss() function. I used the fitDist() on my data and it recommended I used a zero inflated poisson. My response variable is 'deg' and is count data but has a lot of zeros.
fitDEG <- fitDist(deg, data=node_dat, k = 2, type = "counts", try.gamlss = TRUE)
> fitDEG
Family: c("ZIP", "Poisson Zero Inflated")
Fitting method: "nlminb"
Call: gamlssML(formula = y, family = DIST[i])
Mu Coefficients:
[1] 0.3803
Sigma Coefficients:
[1] 2.81
Degrees of Freedom for the fit: 2 Residual Deg. of Freedom 82208
Global Deviance: 38484.9
AIC: 38488.9
SBC: 38507.6
I've tried running a model with a single smoothed term, one numerical explanatory variable (TL), four categorical explanatory variables and two random effects.
mDEG_zip <- gamlss(formula = deg ~ pb(SE_score) + TL + species + sex + season + year +
re(random = ~1|code)+ re(random = ~1|station),
family=ZIP(), data=node_dat)
but I get a warning after twenty iterations
Warning message:
In RS() : Algorithm RS has not yet converged
However I can create a summary output
> summary(mDEG_zip)
******************************************************************
Family: c("ZIP", "Poisson Zero Inflated")
Call: gamlss(formula = deg ~ pb(SE_score) + TL + species + sex + season + year + re(random = ~1 | code) +
re(random = ~1 | station), family = ZIP(), data = node_dat, start.from = mDEG_zip, iter = 20, n.cyc = 40)
Fitting method: RS()
------------------------------------------------------------------
Mu link function: log
Mu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.1461720 0.0755180 -41.661 < 2e-16 ***
pb(SE_score) -0.5060934 0.1431689 -3.535 0.000408 ***
TL 0.0037801 0.0005586 6.767 1.32e-11 ***
speciesSilvertip Shark 2.6530209 0.0326096 81.357 < 2e-16 ***
sexM 0.1816634 0.0277136 6.555 5.60e-11 ***
seasonwet.season -0.0020792 0.0271809 -0.076 0.939026
year2015 0.0614232 0.0449014 1.368 0.171330
year2016 0.1322559 0.0390032 3.391 0.000697 ***
year2017 0.0816437 0.0397759 2.053 0.040115 *
year2018 -0.3669929 0.0557062 -6.588 4.48e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------------------------------
Sigma link function: logit
Sigma Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.28396 0.02217 57.91 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas:
i) Std. Error for smoothers are for the linear effect only.
ii) Std. Error for the linear terms maybe are not accurate.
------------------------------------------------------------------
No. of observations in the fit: 82210
Degrees of Freedom for the fit: 171.7671
Residual Deg. of Freedom: 82038.23
at cycle: 40
Global Deviance: 30763.12
AIC: 31106.66
SBC: 32707.02
******************************************************************
I've tried to use the refit() function but I get the same result after another twenty iterations.
If the model doesn't converge, is this an issue when interpreting the model outputs? A reproducible example dataset is below.
> dput(head(node_dat))
structure(list(station = structure(c(17L, 49L, 23L, 25L, 25L,
9L), .Label = c("BE01", "BE02", "BEUWM01", "BL01", "BL02", "GCB01",
"GCB02", "GCB03", "NI01", "NI01b", "NI03", "PB01", "PB02", "PB03",
"PB04", "PB05", "PB06", "PB07", "PB08", "PB09", "PB10", "PB11",
"PB12", "PB13", "PB14", "PB15", "PB16", "PB17", "PB18", "PB19",
"PB20", "PB21", "PB22", "PB23", "PB24", "PB25", "PB26", "PB27",
"PB28", "PB29", "PB30", "PB4G01", "PB4G02", "PBUWM01", "PBUWM02",
"SA01", "SA02", "SA02b", "SA03", "SA04", "SA05", "SA06", "SA07",
"SA11", "SAUWM01", "SB01", "SB02/AR02", "SB03/AR05", "SB04/AR06",
"VB01", "VB02", "VB03", "VB04"), class = "factor"), monthyear = structure(c(27L,
17L, 38L, 4L, 19L, 29L), .Label = c("2014/01", "2014/02", "2014/03",
"2014/04", "2014/05", "2014/06", "2014/07", "2014/08", "2014/09",
"2014/10", "2014/11", "2014/12", "2015/01", "2015/02", "2015/03",
"2015/04", "2015/05", "2015/06", "2015/07", "2015/08", "2015/09",
"2015/10", "2015/11", "2015/12", "2016/01", "2016/02", "2016/03",
"2016/04", "2016/05", "2016/06", "2016/07", "2016/08", "2016/09",
"2016/10", "2016/11", "2016/12", "2017/01", "2017/02", "2017/03",
"2017/04", "2017/05", "2017/06", "2017/07", "2017/08", "2017/09",
"2017/10", "2017/11", "2017/12", "2018/01", "2018/02", "2018/03",
"2018/04", "2018/05", "2018/06", "2018/07", "2018/08", "2018/09",
"2018/10", "2018/11", "2018/12"), class = "factor"), code = structure(c(99L,
204L, 183L, 146L, 4L, 135L), .Label = c("2390", "13573", "13574",
"13575", "13576", "19318", "19319", "19321", "19322", "19506",
"19514", "19519", "19520", "19524", "25537", "25540", "25541",
"25543", "25546", "25549", "25552", "25553", "27583", "27585",
"27586", "27591", "27592", "27593", "27594", "27595", "27596",
"27597", "27600", "27601", "27605", "27607", "27608", "27613",
"27614", "27617", "27619", "27620", "27621", "27626", "27627",
"27629", "27630", "27631", "27632", "28608", "28611", "28612",
"28618", "28625", "28628", "28629", "28631", "28632", "28633",
"28638", "28641", "28644", "28662", "28672", "28674", "52978",
"54815", "54846", "54852", "54860", "54863", "54865", "54866",
"54868", "54877", "54882", "54883", "54884", "54886", "54890",
"54892", "54895", "54896", "54901", "54904", "54914", "54919",
"54920", "54922", "54925", "54931", "54932", "54938", "54952",
"54954", "54955", "54958", "54959", "54962", "59950", "59953",
"59954", "59955", "59957", "59958", "59959", "59960", "59961",
"59962", "59964", "59966", "59969", "59970", "59971", "59972",
"59973", "59975", "59976", "59979", "59981", "59988", "2388",
"12950", "12952", "12956", "12958", "12960", "12962", "12964",
"12966", "12968", "13577", "14203", "19320", "19523", "25534",
"25535", "25536", "25539", "25542", "25544", "25545", "25547",
"25548", "25550", "27584", "27588", "27589", "27590", "27598",
"27599", "27602", "27603", "27604", "27606", "27609", "27610",
"27611", "27615", "27616", "27618", "27622", "27624", "27625",
"28624", "28627", "28637", "28639", "28642", "28660", "28670",
"34176", "34177", "34178", "34179", "52975", "52977", "54817",
"54821", "54822", "54825", "54845", "54849", "54880", "54887",
"54889", "54893", "54898", "54899", "54905", "54911", "54912",
"54915", "54933", "54947", "54957", "54961", "59951", "59963",
"59968", "59978", "59991", "59992", "59993", "59994", "59995"
), class = "factor"), species = structure(c(1L, 2L, 2L, 2L, 1L,
2L), .Label = c("Grey Reef Shark", "Silvertip Shark"), class = "factor"),
deg = c(0, 0, 0, 0, 0, 0), gs = c(0, 0, 0, 0, 0, 0), btw = c(0,
0, 0, 0, 0, 0), ud = c(0, 0, 0, 0, 0, 0), ri = c(0, 0, 0,
0, 0, 0), SE_score = c(0.35, 0.39, 0.18, 0.23, 0.36, 0.42
), region = structure(c(5L, 6L, 5L, 5L, 5L, 4L), .Label = c("Benares",
"Blenheim", "Grand Chagos Bank", "Nelsons Island", "Peros Banhos",
"Saloman", "Speakers Bank", "Victory Bank"), class = "factor"),
date = structure(c(1456790400, 1430434800, 1485907200, 1396306800,
1435705200, 1462057200), class = c("POSIXct", "POSIXt"), tzone = ""),
month = c(3, 5, 2, 4, 7, 5), season = structure(c(2L, 1L,
2L, 1L, 1L, 1L), .Label = c("dry.season", "wet.season"), class = "factor"),
year = structure(c(3L, 2L, 4L, 1L, 2L, 3L), .Label = c("2014",
"2015", "2016", "2017", "2018"), class = "factor"), sex = structure(c(1L,
2L, 2L, 1L, 1L, 2L), .Label = c("F", "M"), class = "factor"),
TL = c(117, 157, 137, 108, 94, 137), TL_stand = c(0.353383458646617,
0.654135338345865, 0.503759398496241, 0.285714285714286,
0.180451127819549, 0.503759398496241), btw_stand = c(0, 0,
0, 0, 0, 0)), na.action = structure(c(`59` = 59L, `91` = 91L,
`119` = 119L, `144` = 144L, `715` = 715L, `754` = 754L, `780` = 780L,
`803` = 803L, `2116` = 2116L, `2452` = 2452L, `2489` = 2489L,
`2504` = 2504L, `2544` = 2544L, `3070` = 3070L, `3092` = 3092L,
`3126` = 3126L, `3151` = 3151L, `4464` = 4464L, `4800` = 4800L,
`4842` = 4842L, `4862` = 4862L, `4893` = 4893L, `6181` = 6181L,
`8221` = 8221L, `10073` = 10073L, `11232` = 11232L, `11603` = 11603L,
`11639` = 11639L, `11663` = 11663L, `11688` = 11688L, `12266` = 12266L,
`12288` = 12288L, `12322` = 12322L, `12347` = 12347L, `13660` = 13660L,
`14023` = 14023L, `14045` = 14045L, `14075` = 14075L, `14104` = 14104L,
`15417` = 15417L, `15780` = 15780L, `15795` = 15795L, `15837` = 15837L,
`15877` = 15877L, `17138` = 17138L, `17164` = 17164L, `17194` = 17194L,
`17219` = 17219L, `18532` = 18532L, `18895` = 18895L, `18917` = 18917L,
`18951` = 18951L, `18976` = 18976L, `20289` = 20289L, `20652` = 20652L,
`20674` = 20674L, `20704` = 20704L, `20729` = 20729L, `22055` = 22055L,
`22409` = 22409L, `22435` = 22435L, `22461` = 22461L, `22490` = 22490L,
`23803` = 23803L, `24166` = 24166L, `24188` = 24188L, `24218` = 24218L,
`24247` = 24247L, `25560` = 25560L, `25919` = 25919L, `25939` = 25939L,
`25976` = 25976L, `25996` = 25996L, `27308` = 27308L, `27330` = 27330L,
`27360` = 27360L, `27385` = 27385L, `28702` = 28702L, `29065` = 29065L,
`29087` = 29087L, `29121` = 29121L, `29146` = 29146L, `30459` = 30459L,
`30822` = 30822L, `30844` = 30844L, `30874` = 30874L, `30903` = 30903L,
`32216` = 32216L, `32579` = 32579L, `32605` = 32605L, `32631` = 32631L,
`32660` = 32660L, `33973` = 33973L, `34336` = 34336L, `34369` = 34369L,
`34397` = 34397L, `34415` = 34415L, `34875` = 34875L, `34901` = 34901L,
`34931` = 34931L, `34956` = 34956L, `36269` = 36269L, `36632` = 36632L,
`36658` = 36658L, `36684` = 36684L, `36709` = 36709L, `38026` = 38026L,
`38389` = 38389L, `38415` = 38415L, `38441` = 38441L, `38466` = 38466L,
`39783` = 39783L, `40146` = 40146L, `40168` = 40168L, `40198` = 40198L,
`40223` = 40223L, `41540` = 41540L, `41914` = 41914L, `41937` = 41937L,
`41960` = 41960L, `41984` = 41984L, `43297` = 43297L, `43633` = 43633L,
`43675` = 43675L, `43695` = 43695L, `43730` = 43730L, `45014` = 45014L,
`45363` = 45363L, `45400` = 45400L, `45415` = 45415L, `45455` = 45455L,
`45954` = 45954L, `45991` = 45991L, `46009` = 46009L, `46048` = 46048L,
`46541` = 46541L, `46576` = 46576L, `46602` = 46602L, `46627` = 46627L,
`46817` = 46817L, `46859` = 46859L, `46879` = 46879L, `46910` = 46910L,
`48198` = 48198L, `48547` = 48547L, `48589` = 48589L, `48609` = 48609L,
`48640` = 48640L, `49928` = 49928L, `50277` = 50277L, `50319` = 50319L,
`50345` = 50345L, `50370` = 50370L, `51658` = 51658L, `52007` = 52007L,
`52048` = 52048L, `52069` = 52069L, `52100` = 52100L, `53388` = 53388L,
`53737` = 53737L, `53778` = 53778L, `53799` = 53799L, `53830` = 53830L,
`55118` = 55118L, `55467` = 55467L, `55508` = 55508L, `55529` = 55529L,
`55560` = 55560L, `56848` = 56848L, `57197` = 57197L, `57238` = 57238L,
`57264` = 57264L, `57295` = 57295L, `58555` = 58555L, `58596` = 58596L,
`58617` = 58617L, `58648` = 58648L, `59936` = 59936L, `60285` = 60285L,
`60322` = 60322L, `60337` = 60337L, `60377` = 60377L, `60875` = 60875L,
`60905` = 60905L, `60931` = 60931L, `60972` = 60972L, `61441` = 61441L,
`61463` = 61463L, `61497` = 61497L, `61522` = 61522L, `62835` = 62835L,
`63197` = 63197L, `63236` = 63236L, `63260` = 63260L, `63276` = 63276L,
`63793` = 63793L, `64180` = 64180L, `64206` = 64206L, `64232` = 64232L,
`64261` = 64261L, `65574` = 65574L, `65937` = 65937L, `65959` = 65959L,
`65993` = 65993L, `66018` = 66018L, `67331` = 67331L, `67694` = 67694L,
`67716` = 67716L, `67746` = 67746L, `67772` = 67772L, `69088` = 69088L,
`69424` = 69424L, `69466` = 69466L, `69486` = 69486L, `69517` = 69517L,
`70805` = 70805L, `72253` = 72253L, `73419` = 73419L, `73760` = 73760L,
`73802` = 73802L, `73828` = 73828L, `73859` = 73859L, `75141` = 75141L,
`76590` = 76590L, `77752` = 77752L, `78906` = 78906L, `79486` = 79486L,
`79523` = 79523L, `79536` = 79536L, `79556` = 79556L, `80122` = 80122L,
`80159` = 80159L, `80174` = 80174L, `80214` = 80214L, `82125` = 82125L
), class = "omit"), row.names = c(31669L, 80335L, 63799L, 59674L,
1051L, 51949L), class = "data.frame")
You can change the method argument, if you want Rigby and Stasinopoulos Algorithm or Cole and Green, or both, and the number of iteractions. Here is somes examples:
You just need to change the argument
Here the model uses 50 iteractions of RS and 500 of CG. You can use only CG too
Try changing the inicial start values of the parameters, might help. Something like that
But I must warn you, work with random effects in gamlss is quite hard, and is usually that the model doesn't congerge at all, no matter what you do.
Hope this helps