I am using function traitglm() in library(mvabund) to conduct a fourth-corner analysis equivalent to a GLM analysis.
I use the following code:
TraitTGLM = traitglm(AbundanceSO, EnvironmentSO, TraitsSO, method="manyglm", family="negative.binomial",composition = T, col.intercepts = T)
When I run TraitTGLM$fourth, I notice that not all of the environmental columns and trait rows are presented, i.e. this applies to categorical data. I assume that this is a GLM issue and not only a specific TraitGLM issue.
I am following an example, that can be found here: http://rpubs.com/dwarton/68823 under heading Multiple sites - the "fourth corner problem". In this example on ants, none of the expressions of different traits and none of the environmental factors are missing in the example. However, when I replicate the example, I have the same problem with categorical data.
This is my TraitTGLM$fourth output table:
for.L for.O la lo ele soKS soPS
disba 0.71070616 0.685757845 0.77616784 0.44377270 -0.273801917 0.0160309156 -0.32739666
diszo 0.05423857 -0.049717668 0.04232755 -0.07271644 0.172564727 -0.1752525381 -0.02645799
fruca -0.03913508 -0.008395801 -0.10166100 -0.12524280 0.074426535 0.0798772617 -0.02045041
frudr -0.19656983 -0.050675426 0.08987725 0.12030174 -0.143591329 0.1858635739 0.02522805
frufo -0.13801981 0.302674330 0.84252526 -0.96825886 0.062959979 -0.4234008882 0.32380440
frule -0.11621809 -0.052756892 -0.48521160 -0.42152312 0.137112189 0.0539778588 -0.02160870
frunu -0.25605149 -0.154229876 -0.08909875 -0.36192810 -0.009765994 0.1721867978 0.22116788
frusy -0.11173499 -0.050672028 -0.05063261 0.01643027 -0.009587756 0.0579839960 -0.07135602
pgt -0.01447752 0.003283760 0.37786743 0.35312457 -0.106385213 0.0009636494 -0.07629961
polb -0.30055435 -0.198174953 -0.59764188 -0.79383934 -0.057353371 0.3043540174 0.14883309
polf -0.08534353 -0.017679336 -0.02106968 -0.15447035 -0.029297745 0.1438625680 0.01861877
poli -0.35121801 -0.210595674 -0.32407205 -0.68408985 -0.235336244 0.2990251648 0.24455372
polw -0.22586240 -0.148508546 -0.24254181 -0.31744353 0.001497815 0.0464572048 -0.02604469
polx -0.09913082 -0.088810263 -0.30873822 -0.29903480 -0.033386877 0.0907425934 0.06796934
hei 0.07403575 0.021522093 -0.31065910 -0.18872436 0.079194573 -0.0351926077 0.12850278
see -0.14464361 -0.045753203 0.44451589 0.34468547 -0.001742632 -0.0366144001 -0.12092828
woo 0.01117763 0.007088976 0.68168824 0.43018105 -0.152649283 -0.1348766220 -0.02185598
There are 7 columns but it should be 9, i.e. 'for.H' and 'soHS' are missing. Also, there are 17 trait expressions shown but it should be 21. So, those environmental factors and traits with more than 1 possible expression, i.e. 'dis' = 6 expressions ('a','b','f','i','w','x'), one of the expressions is missing.
Could the issue with my data be related to intercepts? I think, the reference level should be '0' and not one of the respective category entries?
Any advice on how to fix this would be much appreciated.
Please find the relevant data for reproducing the analysis below:
AbundanceSO:
EnvironmentSO:
> dput(EnvironmentSO)
structure(list(for. = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L), .Label = c("H", "L", "O"), class = "factor"),
la = c(43.54247057, 43.5453036, 43.54513539, 43.54554091,
43.54569291, 43.54531157, 43.54530949, 43.54578806, 43.54565608,
43.54618175, 43.53344226, 43.53254104, 43.53424692, 43.54350591,
43.53789473, 43.53323841, 43.5386833, 43.53226745, 43.54016387,
43.53381777, 43.53272343, 43.53181684, 43.53344226, 43.53592062,
43.53398407, 43.53394115, 43.53426301, 43.53376949, 43.5280242,
43.52869475, 43.52676892, 43.52626467, 43.52498794, 43.52710688,
43.52534199, 43.52731073, 43.52534199, 43.52087879, 43.52091634,
43.52123821, 43.52058375, 43.52073932, 43.52091634, 43.52135086,
43.52137768, 43.51324522, 43.5133847, 43.51285362, 43.51310575,
43.54387069, 43.54522789, 43.54504013, 43.54276562, 43.54308212,
43.54339325, 43.54300165, 43.54345763, 43.53729928, 43.53766942,
43.53825414, 43.53773379, 43.54011023, 43.5404321, 43.53999758,
43.53366756, 43.53403771, 43.52558875, 43.52663481, 43.52792227,
43.52819049, 43.52402234, 43.52828169, 43.5278579, 43.52444077,
43.52776134, 43.52562094, 43.52537417, 43.52614665, 43.52605546,
43.52589452, 43.52646852, 43.5265168, 43.526034, 43.52467144,
43.52430129, 43.52470362, 43.52629685, 43.52902734, 43.52505231,
43.52439249, 43.52471972, 43.51933384, 43.52024043, 43.52027261,
43.54300165, 43.54258859, 43.54426229, 43.54502404, 43.53538418,
43.54441004, 43.54473972, 43.54309285, 43.54456954, 43.53641415,
43.54307675, 43.54337716, 43.54282999, 43.53600109, 43.5367682,
43.53844726, 43.53772306, 43.53476727, 43.53788936, 43.53802347,
43.53584552, 43.5390588, 43.53507304, 43.53544855, 43.53504622,
43.53284144, 43.53368366, 43.53097999, 43.52645242, 43.53316331,
43.532262, 43.53343153, 43.52754676, 43.52731609, 43.52477336,
43.52390432, 43.52530444, 43.52433348, 43.51094925, 43.51146415,
43.51183438, 43.51266587, 43.51284826, 43.51244056, 43.51255858
), lo = c(122.9251818, 122.9266832, 122.9271353, 122.9272055,
122.9273281, 122.9274899, 122.9276568, 122.9276631, 122.9281064,
122.9286789, 122.928615, 122.9291568, 122.9292373, 122.9296074,
122.9295484, 122.9297844, 122.9298917, 122.9299668, 122.9302565,
122.9301814, 122.9305837, 122.9306696, 122.9308895, 122.9311631,
122.9314957, 122.9334859, 122.933502, 122.9341832, 122.9361144,
122.9365436, 122.9373375, 122.9374019, 122.9379652, 122.9384319,
122.9394726, 122.9399929, 122.9411302, 122.9490212, 122.9492304,
122.9492519, 122.9493806, 122.9495148, 122.9496381, 122.9498527,
122.95077, 122.9510329, 122.9514352, 122.951505, 122.9519717,
122.9225156, 122.9248492, 122.9251603, 122.9262707, 122.9268662,
122.9269627, 122.9273543, 122.9281268, 122.9285184, 122.9288671,
122.9289529, 122.9290173, 122.9293821, 122.9298112, 122.9299454,
122.9335878, 122.9340438, 122.9342744, 122.9343871, 122.9351596,
122.9351864, 122.935181, 122.9356316, 122.9357336, 122.9356906,
122.9357926, 122.935814, 122.9359911, 122.9364685, 122.9365275,
122.9368011, 122.9369298, 122.9371551, 122.9373161, 122.9373429,
122.9377828, 122.937874, 122.9379222, 122.9382066, 122.9384104,
122.938405, 122.9389361, 122.9465697, 122.9479483, 122.9487584,
122.9189161, 122.9236583, 122.9239479, 122.9243127, 122.9243127,
122.924474, 122.9249457, 122.925804, 122.9261636, 122.9261634,
122.9264048, 122.9266141, 122.9266838, 122.9272202, 122.9277138,
122.9278532, 122.9279391, 122.9280034, 122.9285184, 122.9285989,
122.9290334, 122.9292426, 122.9292587, 122.9303155, 122.9311631,
122.9322091, 122.9323272, 122.9330031, 122.9329441, 122.9336522,
122.933974, 122.9340974, 122.9348109, 122.9354332, 122.93627,
122.9367421, 122.9370049, 122.9382924, 122.9469667, 122.9486618,
122.9501853, 122.9504911, 122.9510168, 122.9510275, 122.951403
), ele = c(487, 476, 474, 474, 468, 470, 466, 464, 464, 464,
494, 499, 304, 474, 480, 497, 479, 300, 480, 304, 498, 496,
306, 494, 308, 488, 484, 487, 301, 496, 300, 301, 311, 497,
304, 490, 494, 481, 475, 474, 480, 475, 470, 468, 465, 465,
463, 461, 460, 303, 488, 487, 480, 480, 484, 480, 479, 491,
485, 483, 483, 484, 481, 483, 486, 484, 334, 345, 311, 304,
333, 303, 306, 334, 303, 343, 319, 315, 315, 309, 308, 305,
305, 318, 315, 314, 499, 489, 307, 308, 306, 300, 486, 481,
330, 497, 495, 490, 346, 490, 489, 483, 485, 315, 479, 484,
478, 487, 305, 495, 301, 486, 488, 485, 489, 490, 488, 490,
491, 493, 494, 303, 337, 494, 494, 486, 340, 313, 330, 344,
317, 311, 313, 300, 494, 477, 460, 470, 469), so = structure(c(1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L,
2L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 1L, 3L, 1L,
1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 1L,
1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L, 1L,
1L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L,
3L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L,
1L, 1L, 1L), .Label = c("HS", "KS", "PS"), class = "factor")), .Names = c("for.",
"la", "lo", "ele", "so"), row.names = c("H01", "H02", "H03",
"H04", "H05", "H06", "H07", "H08", "H09", "H10", "H11", "H12",
"H13", "H14", "H15", "H16", "H17", "H18", "H19", "H20", "H21",
"H22", "H23", "H24", "H25", "H26", "H27", "H28", "H29", "H30",
"H31", "H32", "H33", "H34", "H35", "H36", "H37", "H38", "H39",
"H40", "H41", "H42", "H43", "H44", "H45", "H46", "H47", "H48",
"H49", "L01", "L02", "L03", "L04", "L05", "L06", "L07", "L08",
"L09", "L10", "L11", "L12", "L13", "L14", "L15", "L16", "L17",
"L18", "L19", "L20", "L21", "L22", "L23", "L24", "L25", "L26",
"L27", "L28", "L29", "L30", "L31", "L32", "L33", "L34", "L35",
"L36", "L37", "L38", "L39", "L40", "L41", "L42", "L43", "L44",
"L45", "O01", "O02", "O03", "O04", "O05", "O06", "O07", "O08",
"O09", "O10", "O11", "O12", "O13", "O14", "O15", "O16", "O17",
"O18", "O19", "O20", "O21", "O22", "O23", "O24", "O25", "O26",
"O27", "O28", "O29", "O30", "O31", "O32", "O33", "O34", "O35",
"O36", "O37", "O38", "O39", "O40", "O41", "O42", "O43", "O44",
"O45"), class = "data.frame")
TraitsSO:
> dput(TraitsSO)
structure(list(dis = structure(c(1L, 1L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 1L, 1L, 3L, 1L, 1L,
1L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 1L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 1L,
3L, 3L, 3L, 1L, 1L, 3L, 1L), .Label = c("an", "ba", "zo"), class = "factor"),
fru = structure(c(5L, 5L, 3L, 2L, 7L, 3L, 3L, 3L, 2L, 1L,
3L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 6L, 6L, 5L, 2L, 5L, 5L, 1L,
1L, 1L, 3L, 3L, 2L, 4L, 3L, 3L, 2L, 3L, 1L, 3L, 5L, 2L, 7L,
1L, 3L, 3L, 1L, 3L, 3L, 5L, 2L, 3L, 5L, 5L, 5L, 5L, 4L, 2L,
1L, 3L, 3L, 3L, 3L, 3L, 5L), .Label = c("be", "ca", "dr",
"fo", "le", "nu", "sy"), class = "factor"), pg = structure(c(2L,
2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L,
2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L,
2L), .Label = c("s", "t"), class = "factor"), pol = structure(c(2L,
2L, 2L, 4L, 2L, 5L, 4L, 4L, 1L, 4L, 4L, 6L, 4L, 2L, 2L, 4L,
1L, 4L, 1L, 1L, 4L, 4L, 4L, 2L, 4L, 4L, 4L, 2L, 2L, 2L, 4L,
4L, 4L, 2L, 4L, 4L, 4L, 1L, 4L, 4L, 4L, 2L, 4L, 4L, 2L, 2L,
4L, 4L, 2L, 4L, 4L, 2L, 3L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L,
4L), .Label = c("a", "b", "f", "i", "w", "x"), class = "factor"),
hei = c(30, 35, 20, 10, 9, 20, 25, 35, 20, 15, 10, 30, 30,
14, 30, 35, 16, 35, 30, 25, 30, 12, 35, 26, 12, 10, 10, 30,
15, 30, 10, 50, 10, 17, 40, 7, 4, 7, 30, 10, 12, 30, 6, 30,
7, 18, 30, 8, 40, 30, 40, 30, 35, 20, 25, 15, 10, 7, 35,
40, 30, 40), see = c(3, 3, 1, 1, 4, 2, 1, 2, 3, 4, 2, 2,
2, 2, 2, 4, 3, 3, 1, 1, 4, 4, 3, 2, 2, 4, 4, 2, 2, 4, 4,
1, 1, 4, 1, 2, 2, 2, 4, 1, 3, 2, 2, 1, 3, 3, 2, 4, 2, 1,
1, 1, 2, 4, 4, 3, 2, 2, 1, 1, 1, 3), woo = c(1.61, 1.64,
1.61, 1.62, 1.64, 1.61, 1.39, 1.42, 1.73, 1.75, 1.51, 1.77,
1.75, 1.6, 1.88, 1.55, 1.61, 1.72, 1.64, 1.72, 1.65, 1.7,
1.85, 1.74, 1.55, 1.76, 1.76, 1.43, 1.67, 1.48, 1.73, 1.88,
1.4, 1.46, 1.48, 1.67, 1.67, 1.88, 1.53, 1.48, 1.74, 1.64,
1.62, 1.51, 1.64, 1.64, 1.7, 1.64, 1.9, 1.85, 1.7, 1.86,
1.72, 1.43, 1.57, 1.74, 1.52, 1.69, 1.75, 1.45, 1.88, 1.68
)), .Names = c("dis", "fru", "pg", "pol", "hei", "see", "woo"
), row.names = c("albleb", "albodo", "antgha", "apovil", "artlak",
"briret", "buclan", "cansub", "carsph", "catspa", "cropoi", "dalcul",
"dallan", "dalnig", "daloli", "dilobo", "dioehr", "diomal", "dipint",
"diptub", "elltom", "erican", "erysuc", "flesoo", "garcow", "garobt",
"garsoo", "gmearb", "greeri", "halcor", "holpub", "irvoli", "lancor",
"lopdup", "mancal", "memedu", "memscu", "milleu", "mitrot", "morcor",
"ochint", "parama", "pavtom", "permem", "phycol", "phyemb", "ptemac",
"rotwit", "schole", "shoobt", "shorox", "shosia", "sinsia", "stegut",
"steneu", "strnux", "symrac", "syz001", "terala", "tercal", "terche",
"xylxyl"), class = "data.frame")
