[R-sig-ME] seeking input lme4::glmer with a gamma family: link = log or identity?
Scott LaPoint
@dl@point @ending from gm@il@com
Tue Jul 24 16:36:46 CEST 2018
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]]
More information about the R-sig-mixed-models
mailing list