[R-sig-ME] seeking input lme4::glmer with a gamma family: link = log or identity?

Discussion:

Scott LaPoint

2018-07-24 14:36:46 UTC

Hello all,

This is my first posting, so please correct me if I didn't follow all of
the rules. Apologies also for a long post...

I'm using lme4::glmer to explore the influence of several variables on the
distance traveled by birds during their migration. I have chosen a GLMM
because I have some repeated measures (some individuals make multiple
migrations during the study) and some of my potential response variables
are non-normally distributed (although I am only presenting one response
variable here for simplicity).

I have two larger questions/issues where input is welcomed:
1. Is a Gamma distribution best for my distance data? If so, which link
function is most appropriate? I explored two link functions: identity and
log. I have concerns and see potential issues with both (see my annotations
in the reproducible example below.

2. If the log link is the best or most appropriate to use, then the
summary(mDist) produces a sd of the random effect = 0 with the bobyqa
optimizer. Switching to Nelder_Mead gives a reasonable sd, but throws a
convergence warning. What’s going on here and which would be the
appropriate way to move forward? I'm assuming (ick...) that it's because
the log link function is inappropriate, but I'm puzzled why the two
different optimizers produce different sd values.

My concern is not whether there is significance in any of the variables,
but my choice of link function obviously affects my interpretation of the
results!

Any and all help is welcomed! I'm teaching myself GLMMs and would be
thrilled to have input!

#Dataframe:
id <- "101146512 101146512 102057059 102057999 102058563 102058990
102069680 102073750 102073750 102075125 102075125 102075918 10481019
10481021 10481023 10481023 10481023 10481029 150128761 151009981 151010107
151012935 151013389 152198213 152494581 152494617 152494756 16849706
16849887 16849887 16849887 16849913 16849913 16849913 16849913 16849926
171331052 171331052 171331205 171331205 171331516 17750149 17750164
17750164 17750164 17750194 17750194 17750194 17750194 17750233 17750280
17750280 17750280 17750280 183004701 183004701 183004704 19195992 19195992
19196055 19196079 19196079 19196092 19196092 19196092 19196202 19196202
19196202 19196453 19196466 19196466 19196466 19196504 19196504 19196504
212421463 212421463 212421463 214073818 24026401 24026401 24026482 24026482
24026482 24026482 240944289 240944289 240944289 240944308 240944308
240944308 257224604 257224604 257224605 262963605 262963608 262963608
262963608 262963608 262963608 262963609 262963609 262963610 262963611
262963612 262963613 276594658 276594658 45205417 45205417 5009278 5009282
5010023 5010023 5010566 5010567 58674007 58675110 58675110 59530545
59530545 59530620 59530620 59530737 59531771 59533985 59533985 59534620
59534620 60978379 60978551 60978658 60978658 60978784 60978784 6861013
6861013"
id <- strsplit(id, " ")
id <- as.factor(unlist(id))
year <- "2015 2016 2016 2016 2017 2017 2016 2016 2017 2016 2017 2016 2015
2014 2013 2014 2015 2013 2017 2017 2017 2017 2017 2017 2017 2017 2017 2014
2014 2015 2016 2014 2015 2016 2017 2014 2016 2017 2016 2017 2016 2014 2014
2015 2016 2014 2015 2016 2017 2014 2014 2015 2016 2017 2012 2013 2012 2015
2016 2015 2015 2016 2015 2016 2017 2015 2016 2017 2015 2015 2016 2017 2015
2016 2017 2013 2014 2015 2017 2014 2015 2014 2015 2016 2017 2015 2016 2017
2015 2016 2017 2012 2013 2013 2013 2013 2014 2015 2016 2017 2013 2014 2014
2017 2017 2017 2016 2017 2016 2017 2012 2012 2012 2013 2012 2012 2017 2016
2017 2016 2017 2016 2017 2016 2016 2016 2017 2016 2017 2016 2016 2016 2017
2016 2017 2013 2014"
year <- strsplit(year, " ")
year <- as.factor(unlist(year))
age <- "A A SA SA SA SA A SA SA A A A SA SA A A A SA A A A SA A A A A SA SA
SA SA SA SA SA SA A SA SA SA SA A SA A A A A SA SA SA SA A SA SA SA SA SA
SA SA A A A A A A A A SA A A SA SA SA A SA SA SA A A A SA SA SA SA A A A SA
SA A SA SA A A A A SA A A A A A A A SA SA A SA SA SA SA SA SA SA SA SA SA
SA A SA A A A A A A A A A A A A A A A A A SA SA"
age <- strsplit(age, " ")
age <- as.factor(unlist(age))
sex <- "male male female male female male female male male female female
male female male male male male male female male female male female male
female male male female female female female female female female female
female male male female female male female female female female male male
male male female male male male male female female male male male male
female female male male male male male male male male male male female
female female female female female male female female female female female
female male male male female female female female female female male male
male male male male female female male male female male male male male male
female male male male male male male female female male male male male male
male female female female female female female male male male male male
male"
sex <- strsplit(sex, " ")
sex <- as.factor(unlist(sex))
distance <- "3063.84 3077.34 3083.48 2903.85 3385.56 2733.14 2532.30
3619.50 2641.76 3811.18 3626.81 3532.71 3729.12 2729.26 3248.75 3764.56
3397.22 3612.98 2314.65 3372.35 2282.25 2528.39 2218.65 2651.02 3048.75
3198.92 3342.30 2308.78 2766.83 2968.76 2987.37 3686.60 2935.58 2420.89
2314.12 2975.70 4353.93 3895.47 3207.48 2490.26 2816.65 2764.26 2546.95
2622.03 2467.63 3272.38 2819.35 2705.02 2575.35 2518.16 6541.20 6541.10
6562.60 6005.08 3175.18 3531.04 2751.19 2514.32 2762.74 2180.43 2852.11
2586.59 2682.96 2230.56 2019.84 3955.19 4015.25 4238.16 1836.63 4118.42
4059.42 3305.86 3040.52 2670.84 2695.07 3441.04 3350.64 3201.92 2690.18
2798.21 2956.70 2340.82 2571.44 2634.40 2669.62 2914.93 5223.85 2622.16
3644.97 3837.83 3730.73 2992.29 3185.07 3474.46 2904.41 2180.35 1725.05
2779.72 2704.79 2920.02 2702.10 2908.45 3161.65 2990.06 2964.83 2344.53
4519.13 4507.83 3464.48 2090.23 4283.90 5561.30 2364.21 2302.45 2262.61
4666.51 3173.38 2313.64 2470.99 3509.85 3168.33 3380.10 2750.41 2989.91
3451.81 2519.03 2678.53 2370.00 2394.06 1531.50 2339.34 2669.88 2263.27
2245.88 2327.00 3243.36 3859.24"
distance <- strsplit(distance, " ")
distance <- as.numeric(unlist(distance))
direct <- "0.908 0.920 0.891 0.879 0.738 0.929 0.951 0.907 0.914 0.907
0.907 0.918 0.898 0.729 0.868 0.805 0.874 0.839 0.952 0.925 0.943 0.927
0.928 0.938 0.929 0.908 0.911 0.875 0.812 0.775 0.854 0.637 0.750 0.842
0.872 0.810 0.875 0.890 0.802 0.938 0.913 0.943 0.931 0.904 0.878 0.845
0.854 0.859 0.980 0.876 0.849 0.858 0.858 0.899 0.964 0.943 0.948 0.910
0.918 0.896 0.859 0.892 0.857 0.854 0.895 0.902 0.897 0.893 0.910 0.854
0.858 0.898 0.803 0.908 0.886 0.904 0.928 0.944 0.939 0.906 0.878 0.920
0.914 0.923 0.902 0.742 0.407 0.807 0.914 0.895 0.920 0.943 0.944 0.904
0.809 0.731 0.871 0.902 0.926 0.832 0.949 0.926 0.895 0.902 0.913 0.913
0.872 0.863 0.598 0.818 0.736 0.839 0.910 0.883 0.812 0.878 0.810 0.906
0.902 0.877 0.928 0.759 0.868 0.925 0.906 0.912 0.906 0.905 0.901 0.927
0.937 0.893 0.918 0.959 0.943 0.755 0.761"
direct <- strsplit(direct, " ")
direct <- as.numeric(unlist(direct))
CSdirect <- "0.4234405616 0.584627508 0.195104602 0.033922710 -1.859964524
0.705513927 1.001014063 0.410013792 0.504036562 0.410013792 0.410013792
0.557763859 0.289127372 -1.980850943 -0.113827358 -0.960032292 -0.033236412
-0.503350264 1.014445887 0.651786630 0.893559468 0.678650279 0.692082103
0.826400346 0.705513927 0.423445616 0.463741089 -0.019804588 -0.866009522
-1.362987023 -0.301872899 -3.216578784 -1.698782632 -0.463054791
-0.060100061 -0.892873171 -0.019804588 0.181672778 -1.000327765 0.826400346
0.490604738 0.893559468 0.732377576 0.369718318 0.020490885 -0.422759318
-0.301872899 -0.234713777 1.390536969 -0.006372763 -0.369032021
-0.248145602 -0.248145602 0.302559197 1.175627780 0.893559468 0.960718590
0.450309265 0.557763859 0.262263724 -0.234713777 0.208536426 -0.261577426
-0.301872899 0.248831899 0.342854670 0.275695548 0.221968251 0.450309265
-0.301872899 -0.248145602 0.289127372 -0.986895941 0.423445616 0.127945480
0.369718318 0.692082103 0.906991293 0.839832171 0.396581967 0.020490885
0.584627508 0.504036562 0.624922981 0.342854670 -1.806237227 -6.305898385
-0.933168644 0.504036562 0.248831899 0.584627508 0.893559468 0.906991293
0.369718318 -0.906304995 -1.953987295 -0.073531885 0.342854670 0.665218454
-0.597373035 0.974150414 0.665218454 0.248831899 0.342854670 0.490604738
0.490604738 -0.060100061 -0.180986480 -3.740419933 -0.785418576
-1.886828173 -0.503350264 0.450309265 0.087650007 -0.866009522 0.020490885
-0.892873171 0.396581967 0.342854670 0.007059061 0.692082103 -1.577896213
-0.113827358 0.651786630 0.396581967 0.477172913 0.396581967 0.383150143
0.329422845 0.678650279 0.812968522 0.221968251 0.557763859 1.108468658
0.893559468 -1.631623510 -1.551032564"
CSdirect <- strsplit(CSdirect, " ")
CSdirect <- as.numeric(unlist(CSdirect))
start <- "68 59 77 67 94 111 69 79 70 70 73 70 105 107 83 89 90 116 72 75
76 75 67 70 77 73 69 113 120 98 94 104 94 85 92 84 72 70 81 75 91 93 76 68
70 106 87 69 73 82 66 59 69 63 87 89 85 65 65 73 70 69 67 64 71 61 66 61 79
78 76 78 67 65 74 88 87 82 99 76 70 87 78 80 77 74 58 61 83 73 74 87 79 61
76 52 89 69 68 73 66 71 98 103 73 76 92 93 105 121 64 57 54 59 67 70 55 63
73 72 70 62 70 83 67 61 76 65 73 77 75 58 73 68 66 91 89"
start <- strsplit(start, " ")
start <- as.numeric(unlist(start))
CSstart <- "-0.65033704 -1.28876505 -0.01190904 -0.72127349 1.19401052
2.39993009 -0.57940060 0.12996385 -0.50846415 -0.50846415 -0.29565482
-0.50846415 1.97431142 2.11618431 0.41370963 0.83932830 0.91026475
2.75461231 -0.36659126 -0.15378193 -0.08284548 -0.15378193 -0.72127349
-0.50846415 -0.01190904 -0.29565482 -0.57940060 2.54180298 3.03835809
1.47775630 1.19401052 1.90337497 1.19401052 0.55558252 1.05213763
0.48464608 -0.36659126 -0.50846415 0.27183674 -0.15378193 0.98120119
1.12307408 -0.08284548 -0.65033704 -0.50846415 2.04524786 0.69745541
-0.57940060 -0.29565482 0.34277319 -0.79220993 -1.28876505 -0.57940060
-1.00501927 0.69745541 0.83932830 0.55558252 -0.86314638 -0.86314638
-0.29565482 -0.50846415 -0.57940060 -0.72127349 -0.93408282 -0.43752771
-1.14689216 -0.79220993 -1.14689216 0.12996385 0.05902741 -0.08284548
0.05902741 -0.72127349 -0.86314638 -0.22471837 0.76839186 0.69745541
0.34277319 1.54869275 -0.08284548 -0.50846415 0.69745541 0.05902741
0.20090030 -0.01190904 -0.22471837 -1.35970149 -1.14689216 0.41370963
-0.29565482 -0.22471837 0.69745541 0.12996385 -1.14689216 -0.08284548
-1.78532016 0.83932830 -0.57940060 -0.65033704 -0.29565482 -0.79220993
-0.43752771 1.47775630 1.83243853 -0.29565482 -0.08284548 1.05213763
1.12307408 1.97431142 3.10929454 -0.93408282 -1.43063794 -1.64344727
-1.28876505 -0.72127349 -0.50846415 -1.57251083 -1.00501927 -0.29565482
-0.36659126 -0.50846415 -1.07595571 -0.50846415 0.41370963 -0.72127349
-1.14689216 -0.08284548 -0.86314638 -0.29565482 -0.01190904 -0.15378193
-1.35970149 -0.29565482 -0.65033704 -0.79220993 0.98120119 0.83932830"
CSstart <- strsplit(CSstart, " ")
CSstart <- as.numeric(unlist(CSstart))
slat <- "48.93600 48.57967 49.07217 46.62800 46.95150 49.18567 48.38617
42.13933 48.05850 40.65500 41.33650 42.48533 46.64300 47.20000 47.32800
47.67200 46.96500 46.87000 47.76017 48.70567 49.47750 50.51267 52.61100
47.78700 46.23717 47.47950 44.71100 48.23317 48.83300 46.74067 45.30517
47.11000 48.33283 48.97667 49.17733 45.32217 46.29067 48.29750 48.79200
49.32917 46.77550 46.90983 49.07300 48.63250 50.19333 49.04033 49.06967
48.37333 47.25467 48.80133 24.31717 23.90467 24.17783 24.34033 47.08720
44.56210 43.97460 50.48517 48.45483 56.23767 48.57283 50.10050 51.41217
52.84567 53.89317 38.39817 37.69600 37.63700 53.05467 41.47567 42.18283
44.02167 47.89700 48.10050 48.76883 46.12733 46.14600 47.46017 47.44533
50.68483 50.14600 49.12883 48.68533 47.27750 48.00733 33.38517 33.43367
33.40667 45.18433 43.59767 43.66150 49.05917 49.24700 49.96900 47.79283
48.10033 48.90667 41.39850 40.91933 42.05250 47.22350 47.04250 46.11183
50.82667 46.91883 49.67550 39.52133 40.21050 46.11000 46.64850 46.96400
35.70300 37.09417 38.00500 38.98317 34.37200 53.00350 50.92967 49.76167
43.71933 44.08267 48.20750 49.06733 47.95433 43.21917 47.05950 47.03233
49.14417 49.07650 57.30167 48.52517 48.78233 49.98250 51.49683 51.20550
48.39100 48.35700"
slat <- strsplit(slat, " ")
slat <- as.numeric(unlist(slat))
CSs.lat <- "0.50577104 0.44349457 0.52956973 0.10239749 0.15893620
0.54940634 0.40967620 -0.68209590 0.35240870 -0.94151507 -0.82240802
-0.62162482 0.10501907 0.20236702 0.22473782 0.28485936 0.16129562
0.14469229 0.30026899 0.46551583 0.60041002 0.78132865 1.14805777
0.30495812 0.03409139 0.25121577 -0.23264024 0.38293610 0.48776953
0.12208904 -0.12879602 0.18663755 0.40035387 0.51287901 0.54794874
-0.12582490 0.04344169 0.39417918 0.48060388 0.57448611 0.12817635
0.15165346 0.52971479 0.45272777 0.72551699 0.52400499 0.52913280
0.40743213 0.21192180 0.48223451 -3.79690867 -3.86900208 -3.82126139
-3.79286096 0.18265275 -0.25866378 -0.36134227 0.77652242 0.42167604
1.78189778 0.44229913 0.70929292 0.93853598 1.18907150 1.37214506
-1.33594554 -1.45866513 -1.46897667 1.22559882 -0.79808502 -0.67449332
-0.35311576 0.32418303 0.35974912 0.47655442 0.01489444 0.01815744
0.24783742 0.24524381 0.81141738 0.71724504 0.53947230 0.46196097
0.21591184 0.34346562 -2.21207708 -2.20360064 -2.20831948 -0.14991546
-0.42721904 -0.41606335 0.52729769 0.56012510 0.68631041 0.30597704
0.35971940 0.50064498 -0.81157216 -0.89531761 0.69727134 0.20647416
0.17484044 0.01218548 0.83620703 0.15322640 0.63501486 -1.13964873
-1.01920118 0.01186565 0.10598032 0.16112085 -1.80698552 -1.56384810
-1.40466061 -1.23370398 -2.03960692 1.21665574 0.85420853 0.65007495
-0.40595629 -0.34245467 0.37844971 0.52872383 0.33420271 -0.49337021
0.17781157 0.17306301 0.54215331 0.53032649 1.96785509 0.43396950
0.47891384 0.68866983 0.95333217 0.90241587 0.41052035 0.40457811"
CSs.lat <- strsplit(CSs.lat, " ")
CSs.lat <- as.numeric(unlist(CSs.lat))
birds <- data.frame(id = id, year = year, age = age, sex = sex,
distance = distance,
direct = direct, CSdirect = CSdirect,
start = start, CSstart = CSstart,
slat = slat, CSs.lat = CSs.lat)

### Variable definitions:
# id = bird individual id
# year = the year during which the migration took place
# age = age class of individual. Note: some birds transitioned from SA
(subadult) to A (adult) during their tracking lifetime
# sex = sex of the bird
# distance = total distance (km) traveled by bird.id during their spring
migration for that year
# direct = "straightness" index, 0 to 1 = random to straightline, a ratio
of straightline distance from migration start and end point to actual
distance traveled (or: straight/actual)
# CSdirect = centered and scaled 'direct', va base::scale()
# start = Julian day (numeric day of the year) when the migration began
# CSstart = centered and scaled 'start', va base::scale()
# slat = latitude for the start of migration
# CSs.lat = centered and scaled 'slat', va base::scale()

# Most of the continuous variables required scaling for the GLMM and so are
included in the birds dataframe along with the original values.

require(lme4)

# The following models includes the variables that I expect a priori to
influence migration distance.
# Note: the non-scaled variables of slat, direct, and start threw rescaling
warnings.
# First with a Gamma family and an identity link function:
IdB.dist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds,
family = Gamma(link = identity), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(IdB.dist)
# Notes/concerns:
# There are no AIC, etc. estimates provided. The warning messages give
some clues, but how do I address this?
# "Warning messages:
# 1: In vcov.merMod(object, use.hessian = use.hessian) :
# variance-covariance matrix computed from finite-difference Hessian is
# not positive definite or contains NA values: falling back to var-cov
estimated from RX
# 2: In vcov.merMod(object, correlation = correlation, sigm = sig) :
# variance-covariance matrix computed from finite-difference Hessian is
# not positive definite or contains NA values: falling back to var-cov
estimated from RX"
# In general, the estimates for the fixed and random effects look
reasonable. E.g., you might expect a straighter path (high
direct value) would decrease the distance traveled.
# Changing the glmerControl optimizer from "bobyqa" to "Nelder_Mead"
produces similar # estimates and throws the same warnings,
but produces Inf values for AIC, etc. (rather than # NaN):
IdNM.dist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds,
family = Gamma(link = identity), nAGQ = 10, control =
glmerControl(optimizer = "Nelder_Mead", optCtrl=list(maxfun=100000)))
summary(IdNM.dist)

# Now, switching to the log link function:
LogB.dist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = s,
family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(LogB.dist) # sd is 0 for id, so changing to Nelder_Mead produces
accurate random factor sd but a lack of convergence issue...
# Notes/concerns:
# AIC, etc. estimates are now provided.
# There are no warning messages :)
# But, the std deviation for the id (random variable) is now = 0. That
seems highly unlikely.
# Changing the glmerControl optimizer from "bobyqa" to "Nelder_Mead"
produces very different estimates
# and does not converge (but does give a random effects std deviation >
0!):
LogNM.dist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year +
age*sex + (1|id), data = s,
family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "Nelder_Mead", optCtrl=list(maxfun=100000)))
summary(LogNM.dist)

Scott LaPoint
Postdoctoral Researcher, Lamont-Doherty Earth Observatory, Columbia
University
Associate Scientist, Max-Planck Institute for Ornithology
skype: scott_lapoint
twitter @sdlapoint
scottlapoint.weebly.com

[[alternative HTML version deleted]]

Paul Johnson

2018-07-24 17:25:07 UTC

Permalink

Hi Scott,

An incomplete answer…

Post by Scott LaPoint
1. Is a Gamma distribution best for my distance data? If so, which link
function is most appropriate? I explored two link functions: identity and
log. I have concerns and see potential issues with both (see my annotations
in the reproducible example below.

I don’t know (I haven’t run your code) but I’ve always somehow managed to avoid gamma regression for strictly positive data by logging the response and fitting a model with normal errors.

Post by Scott LaPoint
2. If the log link is the best or most appropriate to use, then the
summary(mDist) produces a sd of the random effect = 0 with the bobyqa
optimizer. Switching to Nelder_Mead gives a reasonable sd, but throws a
convergence warning.

(For clarity, I assume that by "sd of the random effect” you mean the square root of the variance parameter that gauges residual inter-bird variation in mean distance and not the SD of the estimate of that parameter, which anyway isn’t output by glmer.)

Why is a random effect variance estimate of zero implausible? I would trust a converged estimate over a non-converged estimate, regardless of whether the estimate is zero. Also… you could compare the log-likelihoods using logLik() — you’d expect the converged fit to have a higher LL. For more general troubleshooting of convergence warnings:
http://rpubs.com/bbolker/lme4trouble1
Another quick check I often do is to fit the non-converged model with glmmTMB (which appears to be more robust than lme4), and compare likelihoods and estimates with lme4.

A quick and dirty model fit assessment is to simulate from the fitted model (which is as easy as simulate(my.fit)), and see if the simulated responses look more or less like the real responses.

Good luck,
Paul

Scott LaPoint

2018-07-25 20:05:44 UTC

Permalink

Thank you Paul, I appreciate your time. And, apologies if my understanding
is often incomplete.

Hi Scott,

Post by Paul Johnson
An incomplete answer…

annotations

Post by Scott LaPoint
in the reproducible example below.

I don’t know (I haven’t run your code) but I’ve always somehow managed to
avoid gamma regression for strictly positive data by logging the response
and fitting a model with normal errors.

If possible, I'd rather not transform the raw data to facilitate
interpretation of the coefficient estimates. I'm likely naive or
misunderstanding something though. Log transforming the distance data does
produce a reasonably normal distribution. The following two models have
very similar AIC, BIC, LogLik, etc. estimates and the p-values of the fixed
effects produce similar interpretations. However, the fixed effects
estimates are quite different.

gammaDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(gammaDist)

logGausDist <- glmer(log(distance) ~ CSs.lat + CSdirect + CSstart + year +
age*sex + (1|id), data = birds, family = gaussian(link = log), nAGQ = 10,
control = glmerControl(optimizer = "bobyqa"))
summary(logGausDist)

The interpretation from these two models are mostly the same: only starting
latitude is a marginally significant predictor of bird migration distance.
Correct?

Post by Paul Johnson
2. If the log link is the best or most appropriate to use, then the

Post by Scott LaPoint
summary(mDist) produces a sd of the random effect = 0 with the bobyqa
optimizer. Switching to Nelder_Mead gives a reasonable sd, but throws a
convergence warning.

Yes, I believe your assumption is correct. In case I am wrong, I'm
referring to these estimates from the summary(model) output:
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.00000 0.0000
Residual 0.02879 0.1697
Number of obs: 137, groups: id, 79

The reason I said that a Std.Dev. = 0 is implausible is because the
ecologist in me says that there is no way that individual birds do not vary
between each other (or even within for birds with multiple migration route
data). Am I misunderstanding the meaning of the Std.Dev here?

Post by Paul Johnson
Another quick check I often do is to fit the non-converged model with
glmmTMB (which appears to be more robust than lme4), and compare
likelihoods and estimates with lme4.
A quick and dirty model fit assessment is to simulate from the fitted
model (which is as easy as simulate(my.fit)), and see if the simulated
responses look more or less like the real responses.
Good luck,
Paul

[[alternative HTML version deleted]]

Thierry Onkelinx

2018-07-25 22:29:15 UTC

Permalink

Dear Scott,

Random effects model only information which is not captured by the fixed
effects. And the random effects are subject to shrinkage. Combine this with
a large number of fixed effect parameters, a small data set and unbalanced
repeated measurements. Then zero variance random effects and convergence
issues don't come as a surprise.

Bottom line: your model is too complex for the data. You'll need to drop
variables or make more observations (often not feasible, I know). Using a
different transformation/link/distribution won't solve any of these issues.

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
***@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>

Post by Scott LaPoint
Thank you Paul, I appreciate your time. And, apologies if my understanding
is often incomplete.
Hi Scott,

Post by Paul Johnson
An incomplete answer…

Post by Scott LaPoint
1. Is a Gamma distribution best for my distance data? If so, which link
function is most appropriate? I explored two link functions: identity

and

Post by Paul Johnson

Post by Scott LaPoint
log. I have concerns and see potential issues with both (see my

annotations

Post by Scott LaPoint
in the reproducible example below.

I don’t know (I haven’t run your code) but I’ve always somehow managed to
avoid gamma regression for strictly positive data by logging the response
and fitting a model with normal errors.

If possible, I'd rather not transform the raw data to facilitate
interpretation of the coefficient estimates. I'm likely naive or
misunderstanding something though. Log transforming the distance data does
produce a reasonably normal distribution. The following two models have
very similar AIC, BIC, LogLik, etc. estimates and the p-values of the fixed
effects produce similar interpretations. However, the fixed effects
estimates are quite different.
gammaDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(gammaDist)
logGausDist <- glmer(log(distance) ~ CSs.lat + CSdirect + CSstart + year +
age*sex + (1|id), data = birds, family = gaussian(link = log), nAGQ = 10,
control = glmerControl(optimizer = "bobyqa"))
summary(logGausDist)
The interpretation from these two models are mostly the same: only starting
latitude is a marginally significant predictor of bird migration distance.
Correct?

Post by Paul Johnson
2. If the log link is the best or most appropriate to use, then the

Post by Scott LaPoint
summary(mDist) produces a sd of the random effect = 0 with the bobyqa
optimizer. Switching to Nelder_Mead gives a reasonable sd, but throws a
convergence warning.

Yes, I believe your assumption is correct. In case I am wrong, I'm
Groups Name Variance Std.Dev.
id (Intercept) 0.00000 0.0000
Residual 0.02879 0.1697
Number of obs: 137, groups: id, 79
The reason I said that a Std.Dev. = 0 is implausible is because the
ecologist in me says that there is no way that individual birds do not vary
between each other (or even within for birds with multiple migration route
data). Am I misunderstanding the meaning of the Std.Dev here?

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

Scott LaPoint

2018-07-26 20:00:40 UTC

Permalink

Thank you Thierry,

The model below does converge and does not produce any warning messages,
but the random effect variance and std dev are both = 0:

mDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex +
(1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))

As I understand it, and please correct me if I'm wrong, it is possible (but
perhaps unlikely) to have these values = 0. If so, I believe this implies
that either the random effect variable is truly not variable or that the
variance of the random effect is being captured by the other fixed effects.
In my case, that might imply that any variation between birds is captured
by year, age, or sex. So, assuming that logic is correct (and it may not
be), then the following model would most likely show a variance and std dev
mDist <- glmer(distance ~ CSs.lat + (1|id), data = birds, family =
Gamma(link = log), nAGQ = 10, control = glmerControl(optimizer = "bobyqa"))

But, it does not, and still shows a variance and std dev of 0. A quick
boxplot of distance grouped by bird id shows both substantial variation
across birds and at times within birds.

Perhaps I'm still missing something? Is 137 observations really too few for
a model with 1 fixed and 1 random effect variable?

Apologies for my ignorance. I do appreciate the guidance while I learn to
swim in the GLMM sea.

scott

Scott LaPoint
Postdoctoral Researcher, Lamont-Doherty Earth Observatory, Columbia
University
Associate Scientist, Max-Planck Institute for Ornithology
skype: scott_lapoint
twitter @sdlapoint
scottlapoint.weebly.com

Post by Thierry Onkelinx
Dear Scott,
Random effects model only information which is not captured by the fixed
effects. And the random effects are subject to shrinkage. Combine this with
a large number of fixed effect parameters, a small data set and unbalanced
repeated measurements. Then zero variance random effects and convergence
issues don't come as a surprise.
Bottom line: your model is too complex for the data. You'll need to drop
variables or make more observations (often not feasible, I know). Using a
different transformation/link/distribution won't solve any of these issues.
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
Havenlaan 88
<https://maps.google.com/?q=Havenlaan+88&entry=gmail&source=g> bus 73,
1000 Brussel
www.inbo.be
////////////////////////////////////////////////////////////
///////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
////////////////////////////////////////////////////////////
///////////////////////////////
<https://www.inbo.be>

Post by Scott LaPoint
Thank you Paul, I appreciate your time. And, apologies if my understanding
is often incomplete.
Hi Scott,

Post by Paul Johnson
An incomplete answer…

Post by Scott LaPoint
1. Is a Gamma distribution best for my distance data? If so, which

link

Post by Paul Johnson

Post by Scott LaPoint
function is most appropriate? I explored two link functions: identity

and

Post by Paul Johnson

Post by Scott LaPoint
log. I have concerns and see potential issues with both (see my

annotations

Post by Scott LaPoint
in the reproducible example below.

I don’t know (I haven’t run your code) but I’ve always somehow managed

Post by Paul Johnson
avoid gamma regression for strictly positive data by logging the

response

Post by Paul Johnson
and fitting a model with normal errors.

If possible, I'd rather not transform the raw data to facilitate
interpretation of the coefficient estimates. I'm likely naive or
misunderstanding something though. Log transforming the distance data does
produce a reasonably normal distribution. The following two models have
very similar AIC, BIC, LogLik, etc. estimates and the p-values of the fixed
effects produce similar interpretations. However, the fixed effects
estimates are quite different.
gammaDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(gammaDist)
logGausDist <- glmer(log(distance) ~ CSs.lat + CSdirect + CSstart + year +
age*sex + (1|id), data = birds, family = gaussian(link = log), nAGQ = 10,
control = glmerControl(optimizer = "bobyqa"))
summary(logGausDist)
The interpretation from these two models are mostly the same: only starting
latitude is a marginally significant predictor of bird migration distance.
Correct?

Post by Paul Johnson
2. If the log link is the best or most appropriate to use, then the

Post by Scott LaPoint
summary(mDist) produces a sd of the random effect = 0 with the bobyqa
optimizer. Switching to Nelder_Mead gives a reasonable sd, but throws

Post by Paul Johnson

Post by Scott LaPoint
convergence warning.

log-likelihoods

Post by Paul Johnson
using logLik() — you’d expect the converged fit to have a higher LL.

For

Post by Paul Johnson
http://rpubs.com/bbolker/lme4trouble1

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

Thierry Onkelinx

2018-07-27 08:19:14 UTC

Permalink

Dear Scott,

If you run the model with Laplace approximation instead of Adaptive
Gauss-Hermite Quadrature, then the random effect yields has sensible
variance estimate. Maybe someone else (Ben Bolker?) can chime in and
explain why.

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
***@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>

Post by Scott LaPoint
Thank you Thierry,
The model below does converge and does not produce any warning messages,
mDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex +
(1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
As I understand it, and please correct me if I'm wrong, it is possible
(but perhaps unlikely) to have these values = 0. If so, I believe this
implies that either the random effect variable is truly not variable or
that the variance of the random effect is being captured by the other fixed
effects. In my case, that might imply that any variation between birds is
captured by year, age, or sex. So, assuming that logic is correct (and it
may not be), then the following model would most likely show a variance and
mDist <- glmer(distance ~ CSs.lat + (1|id), data = birds, family =
Gamma(link = log), nAGQ = 10, control = glmerControl(optimizer = "bobyqa"))
But, it does not, and still shows a variance and std dev of 0. A quick
boxplot of distance grouped by bird id shows both substantial variation
across birds and at times within birds.
Perhaps I'm still missing something? Is 137 observations really too few
for a model with 1 fixed and 1 random effect variable?
Apologies for my ignorance. I do appreciate the guidance while I learn to
swim in the GLMM sea.
scott
Scott LaPoint
Postdoctoral Researcher, Lamont-Doherty Earth Observatory, Columbia
University
Associate Scientist, Max-Planck Institute for Ornithology
skype: scott_lapoint
scottlapoint.weebly.com
On Wed, Jul 25, 2018 at 6:29 PM, Thierry Onkelinx <

Post by Thierry Onkelinx
Dear Scott,
Random effects model only information which is not captured by the fixed
effects. And the random effects are subject to shrinkage. Combine this with
a large number of fixed effect parameters, a small data set and unbalanced
repeated measurements. Then zero variance random effects and convergence
issues don't come as a surprise.
Bottom line: your model is too complex for the data. You'll need to
drop variables or make more observations (often not feasible, I know).
Using a different transformation/link/distribution won't solve any of
these issues.
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
Havenlaan 88
<https://maps.google.com/?q=Havenlaan+88&entry=gmail&source=g> bus 73,
1000 Brussel
www.inbo.be
////////////////////////////////////////////////////////////
///////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
////////////////////////////////////////////////////////////
///////////////////////////////
<https://www.inbo.be>

Post by Scott LaPoint
Thank you Paul, I appreciate your time. And, apologies if my
understanding
is often incomplete.
Hi Scott,

Post by Paul Johnson
An incomplete answer…

Post by Scott LaPoint
1. Is a Gamma distribution best for my distance data? If so, which

link
identity and

Post by Paul Johnson

Post by Scott LaPoint
log. I have concerns and see potential issues with both (see my

annotations

Post by Scott LaPoint
in the reproducible example below.

I don’t know (I haven’t run your code) but I’ve always somehow managed

Post by Paul Johnson
avoid gamma regression for strictly positive data by logging the

response

Post by Paul Johnson
and fitting a model with normal errors.

If possible, I'd rather not transform the raw data to facilitate
interpretation of the coefficient estimates. I'm likely naive or
misunderstanding something though. Log transforming the distance data does
produce a reasonably normal distribution. The following two models have
very similar AIC, BIC, LogLik, etc. estimates and the p-values of the fixed
effects produce similar interpretations. However, the fixed effects
estimates are quite different.
gammaDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(gammaDist)
logGausDist <- glmer(log(distance) ~ CSs.lat + CSdirect + CSstart + year +
age*sex + (1|id), data = birds, family = gaussian(link = log), nAGQ = 10,
control = glmerControl(optimizer = "bobyqa"))
summary(logGausDist)
The interpretation from these two models are mostly the same: only starting
latitude is a marginally significant predictor of bird migration distance.
Correct?

Post by Paul Johnson
2. If the log link is the best or most appropriate to use, then the

Post by Scott LaPoint
summary(mDist) produces a sd of the random effect = 0 with the bobyqa
optimizer. Switching to Nelder_Mead gives a reasonable sd, but

throws a

Post by Paul Johnson

Post by Scott LaPoint
convergence warning.

log-likelihoods

Post by Paul Johnson
using logLik() — you’d expect the converged fit to have a higher LL.

For

Post by Paul Johnson
http://rpubs.com/bbolker/lme4trouble1

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]

Scott LaPoint

2018-07-27 18:14:09 UTC

Permalink

Thank you again Thierry. I think I had explored different nAGQ values
previously, but not very systematically.

Following Thierry's idea and Professor Bolker's suggestion elsewhere (here:
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html) under "Trouble
shooting", I compared the LogLikelihood estimates for models run with
different model fitting options:

nAGQ = 0 ("optimizing the random effects and the fixed-effects coefficients
in the penalized iteratively reweighted least squares step")
nAGQ = 1 (Laplace)
nAGQ = 10 (>1 = Gauss-Hermite)

I run this model, with three variants of the nAGQ argument:
mDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex +
(1|id), data = birds, family = Gamma(link = log), nAGQ = ?, control =
glmerControl(optimizer = "bobyqa"))

I used Bolker's suggestion:
mods <- allFit(mDist)

I provide AIC values as well, just for reference.

With nAGQ = 0:
AIC = 1887.4
allFit() = 100% of the optimizers produce models that converge. The LogLik
values all match at -929.6804

With nAGQ = 1:
AIC = 1843.5
allFit() = 50% of the optimizers produce models that converge. The LogLik
values for each model are very similar (range = -908.2571 to -907.7478)

With nAGQ = 10:
AIC = 31.8
allFit() = 50% of the optimizers produce models that converge. The LogLik
values vary much more substantially (range = -10.481851 to -1.879633)

So, I'm puzzled why despite increasing the nAGQ argument, I seem to be
achieving less certainty...

Thanks to you all.

Scott LaPoint
Postdoctoral Researcher, Lamont-Doherty Earth Observatory, Columbia
University
Associate Scientist, Max-Planck Institute for Ornithology
skype: scott_lapoint
twitter @sdlapoint
scottlapoint.weebly.com

Post by Thierry Onkelinx
Dear Scott,
If you run the model with Laplace approximation instead of Adaptive
Gauss-Hermite Quadrature, then the random effect yields has sensible
variance estimate. Maybe someone else (Ben Bolker?) can chime in and
explain why.
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
Havenlaan 88
<https://maps.google.com/?q=Havenlaan+88&entry=gmail&source=g> bus 73,
1000 Brussel
www.inbo.be
////////////////////////////////////////////////////////////
///////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
////////////////////////////////////////////////////////////
///////////////////////////////
<https://www.inbo.be>

Post by Thierry Onkelinx
Dear Scott,
Random effects model only information which is not captured by the fixed
effects. And the random effects are subject to shrinkage. Combine this with
a large number of fixed effect parameters, a small data set and unbalanced
repeated measurements. Then zero variance random effects and convergence
issues don't come as a surprise.
Bottom line: your model is too complex for the data. You'll need to
drop variables or make more observations (often not feasible, I know).
Using a different transformation/link/distribution won't solve any of
these issues.
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
Havenlaan 88
<https://maps.google.com/?q=Havenlaan+88&entry=gmail&source=g> bus 73,
1000 Brussel
www.inbo.be
////////////////////////////////////////////////////////////
///////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
////////////////////////////////////////////////////////////
///////////////////////////////
<https://www.inbo.be>

Post by Scott LaPoint
Thank you Paul, I appreciate your time. And, apologies if my understanding
is often incomplete.
Hi Scott,

Post by Paul Johnson
An incomplete answer…

Post by Scott LaPoint
1. Is a Gamma distribution best for my distance data? If so, which

link
identity and

Post by Paul Johnson

Post by Scott LaPoint
log. I have concerns and see potential issues with both (see my

annotations

Post by Scott LaPoint
in the reproducible example below.

I don’t know (I haven’t run your code) but I’ve always somehow

managed to

Post by Paul Johnson
avoid gamma regression for strictly positive data by logging the

response

Post by Paul Johnson
and fitting a model with normal errors.

If possible, I'd rather not transform the raw data to facilitate
interpretation of the coefficient estimates. I'm likely naive or
misunderstanding something though. Log transforming the distance data does
produce a reasonably normal distribution. The following two models have
very similar AIC, BIC, LogLik, etc. estimates and the p-values of the fixed
effects produce similar interpretations. However, the fixed effects
estimates are quite different.
gammaDist <- glmer(distance ~ CSs.lat + CSdirect + CSstart + year + age*sex
+ (1|id), data = birds, family = Gamma(link = log), nAGQ = 10, control =
glmerControl(optimizer = "bobyqa"))
summary(gammaDist)
logGausDist <- glmer(log(distance) ~ CSs.lat + CSdirect + CSstart + year +
age*sex + (1|id), data = birds, family = gaussian(link = log), nAGQ = 10,
control = glmerControl(optimizer = "bobyqa"))
summary(logGausDist)
The interpretation from these two models are mostly the same: only starting
latitude is a marginally significant predictor of bird migration distance.
Correct?

Post by Paul Johnson
2. If the log link is the best or most appropriate to use, then the

Post by Scott LaPoint
summary(mDist) produces a sd of the random effect = 0 with the

bobyqa

Post by Paul Johnson

Post by Scott LaPoint
optimizer. Switching to Nelder_Mead gives a reasonable sd, but

throws a

Post by Paul Johnson

Post by Scott LaPoint
convergence warning.

Post by Paul Johnson
whether the estimate is zero. Also… you could compare the

log-likelihoods

Post by Paul Johnson
using logLik() — you’d expect the converged fit to have a higher LL.

For

Post by Paul Johnson
http://rpubs.com/bbolker/lme4trouble1

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

[[alternative HTML version deleted]]