I have a data set where two traits (lambda and lFC.inf) are measured simultaneously in different species, and the same species may be observed in more than one experiment. I apply a bivariate model to these data, where I want to check whether the two traits are correlated across species (e.g. does species that have a large value for one trait have a small value for the other). For this, I compare two models with struct ="UN" or struct = "DIAG", using AICc. This gives strong support for the model with struct = "UN", with an estimated negative correlation of -0.71.
I want to visualize this pattern, and naively thought that I could just obtain the species-level random effects for lambda and lFC.inf from this model (obtained using ranef) and make a plot between them. However, this plot shows a very (unrealistic) tight relationship, and the negative correlation between these is considerably stronger (-0.94) than what the model estimated. Interestingly, when I do the same for the model with struct = "DIAG", the correlation becomes more similar to the one estimated in the model (-0.67).
So my questions are (1) is my interpretation of the model comparison correct? and (2) how do I best visualize the result in a figure?
See below for script.
Any input on this would be greatly appreciated.
Best,
Sigurd Einum
library(dplyr)
library(tidyr)
library(MuMIn)
library(metafor)
data <- structure(list(species = c("Acanthopagrus schlegelii", "Acanthopagrus schlegelii",
"Acipenser persicus", "Acipenser persicus", "Acipenser sinensis",
"Acipenser transmontanus", "Alosa alosa", "Anguilla japonica",
"Anguilla japonica", "Anguilla marmorata", "Anguilla marmorata",
"Anguilla marmorata", "Anguilla marmorata", "Apeltes quadracus",
"Callinectes ornatus", "Callinectes sapidus", "Callinectes sapidus",
"Carcinus maenas", "Carcinus maenas", "Carcinus maenas", "Carcinus maenas",
"Chasmagnathus granulata", "Chasmagnathus granulata", "Dichotomyctere nigroviridis",
"Eleginops maclovinus", "Eleginops maclovinus", "Eleginops maclovinus",
"Eleginops maclovinus", "Eleginops maclovinus", "Eleginops maclovinus",
"Eleginops maclovinus", "Eriocheir sinensis", "Eriocheir sinensis",
"Fundulus heteroclitus", "Fundulus heteroclitus", "Fundulus heteroclitus",
"Galaxias maculatus", "Gasterosteus aculeatus", "Gasterosteus aculeatus",
"Hemigrapsus nudus", "Heteropneustes fossilis", "Leptograpsus variegatus",
"Macrobrachium acanthurus", "Macrobrachium amazonicum", "Morone saxatilis",
"Morone saxatilis", "Morone saxatilis", "Morone saxatilis", "Oncorhynchus keta",
"Oncorhynchus mykiss", "Oncorhynchus mykiss", "Oncorhynchus mykiss",
"Oreochromis mossambicus", "Oreochromis mossambicus", "Oreochromis mossambicus",
"Oreochromis mossambicus", "Oreochromis mossambicus", "Oreochromis mossambicus",
"Oreochromis niloticus", "Palaemon northropi", "Paralichthys lethostigma",
"Paralichthys lethostigma", "Penaeus vannamei", "Penaeus vannamei",
"Penaeus vannamei", "Portunus trituberculatus", "Portunus trituberculatus",
"Portunus trituberculatus", "Salmo salar", "Salmo salar", "Salmo salar",
"Salmo salar", "Salmo trutta", "Salmo trutta", "Salmo trutta",
"Salmo trutta", "Salmo trutta", "Scylla paramamosain", "Scylla paramamosain",
"Solea senegalensis", "Solea senegalensis", "Solea senegalensis",
"Solea senegalensis", "Solea senegalensis", "Sparus auratus",
"Sparus auratus", "Sparus auratus", "Synechogobius ommaturus",
"Synechogobius ommaturus", "Synechogobius ommaturus", "Synechogobius ommaturus",
"Trachidermus fasciatus", "Uca pugilator", "Uca pugnax"), lambda = c(1.05112521267463,
0.271868325477894, 0.0298790410821654, 0.00534951064664712, 0.0711554418256704,
0.019557352594906, 0.0165994917942084, 0.00681365758525864, 0.00610335877546237,
0.029899996039646, 0.143755757058074, 3.00332717555465, 2.6730872261002,
0.00612842450151681, 0.00109519468407754, 0.00047742265820241,
0.00322281473823271, 0.104074230037889, 0.0288707198573549, 0.00354404467661802,
0.000749756241195731, 0.00215762492662288, 0.00499787464672858,
0.165201377203689, 0.00684993821105755, 0.00981038995299565,
0.0171242626131419, 0.00579552816584897, 0.0262022925832092,
0.0873615180472193, 0.00192902545011328, 0.0660321771320779,
0.0672058221676321, 0.0130762009465099, 0.356466412662354, 0.0151721601733623,
0.13018440831611, 0.00387222800690919, 0.00774666799294241, 0.0063393212586308,
0.0025769270054248, 0.00659147689579521, 0.00543906446322283,
3.78253200529015, 0.00920754678220984, 0.0180655524505748, 0.0100895104374508,
0.00974361944957595, 0.00511956696512464, 0.347064925263177,
0.0519871367708311, 0.473255575900437, 0.0942846149603939, 0.00165924542650299,
0.0424763415887829, 0.150727301167667, 2.04742920134463, 0.0665216185603376,
0.0117666739554029, 0.697695078462888, 0.0205833447075628, 0.0143042522170995,
0.260866287723015, 1.07007078109886, 0.176736704265223, 0.0549845026250279,
0.0342134266003462, 0.0336203779355841, 0.00560394242453064,
0.139913260837712, 0.0205997635334676, 0.0087003349127729, 0.000186242562636179,
0.00381649487838412, 0.00876175685041844, 0.00650242635796131,
0.00824503845878259, 0.0702320548363362, 0.0154919642853277,
0.00543746352013002, 0.0033432113872422, 0.00258097906789658,
0.00490541064577789, 0.00495077159457755, 0.0119285569104972,
0.00784011551815093, 0.00896171207032962, 6.77872590968513, 0.141160214515369,
0.177449533069521, 0.0365211378186829, 0.0515516999629714, 0.0100223576547731,
0.0132738582293963), lambda.var = c(0.229770805296645, 0.00311319484562275,
8.30335137134461e-05, 1.39833144739391e-05, 0.000511409666073228,
3.35698268886982e-06, 2.34436807086907e-05, 3.9461516157468e-06,
6.17777694171228e-06, 0.000225162839047296, 0.000711042638047283,
6.86463335142108, 33.7993683348458, 8.69467655862716e-06, 5.34673414495153e-05,
8.84385829663742e-07, 6.45509984032976e-06, 0.0504639381166314,
0.000264137620323017, 1.34549353214643e-06, 8.79669188905294e-07,
5.50079444210757e-06, 2.08983299415301e-06, 0.0100899078144925,
2.31764867431424e-06, 7.25002515095353e-05, 3.26335098888596e-05,
3.36750464855695e-05, 1.42667906004078e-05, 0.00586284841885338,
8.4006692969467e-05, 0.00107299646850373, 0.00108464579223292,
1.37872452418955e-05, 0.0254429269051615, 1.5554145410443e-05,
0.000211148554291663, 4.89868407769818e-07, 4.58538383930258e-05,
1.43093743310837e-05, 1.69975345633274e-05, 1.73221142336391e-05,
5.27439349253101e-05, 888.055253438518, 5.14529301990942e-05,
0.000165322556459889, 5.22435865898896e-05, 7.92005116060282e-06,
3.49839411947353e-06, 0.22933072472968, 0.00149753490550964,
1.17282138460607, 0.00134169772237486, 0.000392170309172174,
0.000436208848455048, 0.00113263578586526, 0.00016310318295376,
0.000249498952286906, 1.43363341989909e-05, 0.470544076803978,
0.000149010597825891, 1.09123575439209e-05, 0.000180061655585789,
0.00743560506046289, 0.000547267889577824, 0.000432087999023795,
0.000587581079958796, 6.99788967524493e-05, 5.40317927547147e-06,
0.0174195268154505, 1.9189557706461e-05, 3.75026849669255e-05,
5.50108353758318e-06, 1.46979827359604e-07, 5.05070782243564e-06,
1.75618457380686e-06, 2.27561742168883e-06, 0.000742744603066814,
7.04618305671455e-05, 2.98610943596692e-06, 4.02035778600495e-06,
1.66642016305517e-06, 1.86653739386804e-06, 6.84780350253536e-06,
7.38532164970507e-07, 5.26205419435704e-05, 2.06420465600515e-06,
19.8682611742585, 0.00204261459750374, 0.00632813919114431, 0.000393082688978805,
2.42631856269185e-06, 1.03934972880783e-05, 2.14971337739398e-06
), lFC.inf = c(0.00865303841411384, 0.0498435498654334, 0.115754250451254,
0.32328798532204, 0.110026621936279, 0.103903022240831, 0.0235083982965177,
0.0850591611407937, 0.0710273910467964, 0.126029408424364, 0.109637302451021,
0.0364460179579933, 0.0323762434313177, 0.0362182591130537, 0.459997826336228,
0.160782563735025, 0.0228369749686444, 0.023002274937268, 0.0641812301605438,
0.0539178731743926, 0.0948944419484686, 0.0672727487932266, 0.0562211463298234,
0.0270178723804437, 0.0585191089320559, 0.0342365832356017, 0.0190707136053501,
0.0387175599579897, 0.0148765533754871, 0.0138588759272592, 0.111345179666044,
0.137255479129647, 0.0985749310605587, 0.0196506549099298, 0.0138758292739683,
0.0341012897232172, 0.0343441691528201, 0.0562499426501846, 0.0308832826270344,
0.0321235090530118, 0.252758106164492, 0.0475269653408294, 0.187845383849484,
0.0503219265909209, 0.0436774856398576, 0.0084041169791645, 0.0607433110135273,
0.0256481428220131, 0.0276988097727251, 0.017069599929371, 0.0280627653917539,
0.0201389816306086, 0.0279970175946409, 0.229244683198107, 0.0842064919081034,
0.064653569316684, 0.0336771555509742, 0.0527868895390138, 0.155981060350519,
0.0848631229565141, 0.051475410360614, 0.0257119866932405, 0.746443890083462,
0.317122891992046, 1.57636928253702, 0.0685164156293337, 0.254244957442793,
0.229936415585668, 0.0473357549099458, 0.00955401476742476, 0.0447964002463933,
0.067662359021822, 0.89582053073337, 0.110625335908012, 0.101801862884989,
0.0289752012368463, 0.111760108801887, 0.0735175467878335, 0.0790681248870101,
0.0547838067415764, 0.0410141844461181, 0.0480599925580839, 0.0534305309956142,
0.0270334269258641, 0.0224425263595611, 0.00847866439878538,
0.0436477040099897, 0.0448743539438881, 0.0302574213074796, 0.010003424768621,
0.0567043312447577, 0.0117732690178457, 0.0177621566361039, 0.0396200858957789
), lFC.inf.var = c(4.76095901086041e-07, 4.63277944313332e-06,
0.000154344126395168, 0.0187513180708446, 7.85792978102002e-05,
7.40710634781114e-06, 4.73908816254375e-06, 0.000149211720342422,
0.000160662464131529, 0.000462911850931987, 2.26348473013204e-05,
2.52202369432042e-06, 1.9050526903333e-05, 7.13282047388166e-05,
7.46010718395484, 0.0854836307557864, 0.000111307767777283, 0.000125863442354,
0.000121696853988268, 8.10198252546276e-05, 0.010261696249353,
0.00165923107006046, 3.88657402243703e-05, 1.94601624021025e-05,
3.44668212584297e-05, 0.000127528495785826, 3.72272462502041e-06,
0.000362272167779435, 3.11214198036298e-07, 2.60546801017905e-06,
0.161526839025964, 0.000512235350973215, 0.000254630341404198,
6.69955938927087e-06, 1.34821804647799e-06, 6.92021550477553e-06,
4.25032856214131e-07, 9.97469539608499e-06, 6.71905011305909e-05,
9.33541582618414e-05, 0.0783052191608477, 0.00014818066740766,
0.0231580781811404, 0.000235204878302195, 0.000290283283062689,
5.36012957203518e-06, 0.00043461381693984, 7.32970148100876e-06,
2.60709046368928e-05, 2.68427105605434e-05, 3.5231295436604e-05,
6.92043684552128e-05, 1.15249124831978e-05, 6.99306220439793,
0.000253059999366038, 7.54713588429824e-06, 5.7416068667544e-10,
1.8327542326458e-05, 0.000705673570495417, 0.000176916217447162,
9.15161780840805e-05, 5.0311495592456e-06, 4.44997043928075e-05,
1.36951971968623e-05, 0.00206688520259364, 8.90067365545625e-05,
0.00639899419045084, 0.000656552884753602, 3.88841877774328e-05,
2.24576033096961e-06, 8.86119375567365e-06, 0.000709675849256686,
120.3957975123, 1.5933146958773e-05, 9.79592585950157e-05, 5.85300144218265e-06,
6.3943462197904e-05, 4.90089718212599e-05, 0.000143560250039281,
6.56474930775405e-05, 0.000203192432402931, 0.000240200429497597,
5.25703038724162e-05, 4.80702019940201e-05, 3.18771246695284e-07,
1.09723519205789e-05, 7.68412733822801e-06, 9.58106036631253e-07,
6.10047499524103e-06, 1.03899573602868e-06, 7.06308100644178e-05,
2.30561078270423e-08, 3.09969381771988e-06, 1.51615556548692e-06
), lFC.lambda.covar = c(-0.000107218275084689, -6.86024367942994e-05,
-7.93067406991034e-05, -0.000495817329492644, -0.00010328031536203,
-3.73611836829387e-06, -7.32820957582603e-06, -2.23104125232692e-05,
-2.31051333424716e-05, -0.000240513336205618, -6.15877995918089e-05,
-0.00133555978593521, -0.00805475627421848, -2.32236087643178e-05,
-0.0199624091807527, -0.000274886078971075, -2.62237942591558e-05,
-0.00127223362678345, -0.000113034567671046, -9.96844146728906e-06,
-9.49029115171746e-05, -9.17170156032077e-05, -7.24885045812932e-06,
-0.000256604606135574, -8.11011751237583e-06, -7.92891830787433e-05,
-7.52328043746124e-06, -0.000103302683734247, -1.26171738235285e-06,
-5.86858635947634e-05, -0.00366663154067903, -0.000538796075191493,
-0.000379555074330538, -8.6821477600903e-06, -7.30514400268722e-05,
-7.13617541580289e-06, -3.99339507664623e-06, -1.54774348180237e-06,
-2.99564308191546e-05, -3.42848522309531e-05, -0.00114309150642882,
-3.98478699648329e-05, -0.00107695195364669, -0.200641831743645,
-0.000114791368203266, -2.47635857286621e-05, -0.000139311354936841,
-5.46544646075489e-06, -8.24434687620744e-06, -0.00125883320678818,
-0.000134158292699473, -0.00466125459986037, -7.98662870680556e-05,
-0.0523657589759199, -0.000236464677972336, -6.36426559039074e-05,
-1.40353147667711e-07, -4.69022936429114e-05, -9.30677738894666e-05,
-0.00306148362921299, -8.78345622981773e-05, -6.45699757049734e-06,
-4.54156013605559e-05, -8.62477867411103e-05, -0.00065605894250463,
-0.000151780624088738, -0.00174424909983081, -0.000193593133627035,
-1.15152712312432e-05, -7.40487586337882e-05, -8.31788745985921e-06,
-0.000157062070468861, -0.0257346930080777, -7.97914207216513e-07,
-1.61480375881298e-05, -2.38834383063823e-06, -6.06851023805994e-06,
-8.97968117979233e-05, -5.93538458904149e-05, -1.26280099602993e-05,
-2.76646141718886e-05, -1.96712747117823e-05, -9.12256931766496e-06,
-1.66798946997085e-05, -3.74920728211624e-07, -2.11326806992407e-05,
-3.3762564879086e-06, -0.00130285017473134, -5.75300750638565e-05,
-3.83062378694734e-05, -8.81682107886266e-05, -2.13546807644931e-07,
-3.9656757057916e-06, -1.13134245574533e-06)), row.names = c("4",
"10", "22", "23", "24", "26", "30", "36", "37", "40", "42", "43",
"44", "45", "47", "57", "65", "66", "67", "77", "78", "80", "83",
"122", "133", "137", "140", "141", "142", "143", "144", "179",
"180", "186", "193", "198", "199", "201", "203", "207", "212",
"214", "215", "218", "219", "220", "222", "223", "227", "234",
"235", "236", "238", "242", "243", "246", "247", "248", "249",
"255", "258", "259", "261", "263", "271", "280", "283", "284",
"285", "287", "291", "293", "298", "299", "301", "302", "304",
"337", "338", "376", "381", "382", "385", "388", "402", "413",
"432", "452", "453", "454", "455", "461", "462", "463"), class = "data.frame")
data$id <- 1:nrow(data)
# Create a new dataframe
new_data <- data.frame(
yi = c(data$lambda, data$lFC.inf),
species = rep(data$species, 2),
id = rep(data$id, 2),
lambda.var = rep(data$lambda.var,2),
lFC.inf.var = rep(data$lFC.inf.var,2),
lFC.lambda.covar =rep(data$lFC.lambda.covar,2)
)
new_data$dep.trait <- rep(c("lambda", "lFC"), each = nrow(data))
new_data$v1i <- ifelse(new_data$dep.trait == "lambda",new_data$lambda.var,new_data$lFC.lambda.covar)
new_data$v2i <- ifelse(new_data$dep.trait == "lFC",new_data$lFC.inf.var,new_data$lFC.lambda.covar)
# Sort the DataFrame by id in ascending order and trait in reverse alphabetical order
new_data <- new_data %>%
arrange(id, desc(dep.trait))
#Create variance-covariance matrix for the errors in the parameter estimates (lambda and lFC.inf)
V <- vcalc(vi=1, cluster=id, rvars=c(v1i, v2i), data=new_data)
#allow for different variances and correlation between the dependent traits,
mod1 <- rma.mv(yi = yi, V = V, mods = ~dep.trait, data = new_data,
random = list(~ dep.trait| species, ~1|id),
struct ="UN", method = "ML")
#allow for different variances for each dependent trait but no correlation,
mod2 <- rma.mv(yi = yi, V = V, mods = ~dep.trait, data = new_data,
random = list(~ dep.trait | species, ~1|id),
struct ="DIAG", method = "ML")
AICc(mod1,mod2)
summary(mod1)
#extract species level random effects from mod1
coef <- ranef(mod1)
species.coef <- as.data.frame(coef$`~dep.trait | species`)
species.coef$names <- rownames(species.coef)
rownames(species.coef) <- NULL
species.coef <- species.coef %>%
separate(names, into = c("parameter", "names"), sep = "\\|")
species.coef$parameter <- trimws(species.coef$parameter)
species.coef$names <- trimws(species.coef$names)
head(species.coef)
species.coef <- species.coef[,c(1,5,6)]
str(species.coef)
species.coef <- species.coef %>%
pivot_wider(names_from = parameter, values_from = intrcpt)
species.coef <- as.data.frame(species.coef)
par(mfrow= c(2,1))
#correlation between lambda and lFC.inf across species based on random effects from mod1
cor.test(species.coef$lambda, species.coef$lFC)
plot(species.coef$lambda, species.coef$lFC, main = "Model 1")
#extract species level random effects from mod2
coef <- ranef(mod2)
species.coef <- as.data.frame(coef$`~dep.trait | species`)
species.coef$names <- rownames(species.coef)
rownames(species.coef) <- NULL
species.coef <- species.coef %>%
separate(names, into = c("parameter", "names"), sep = "\\|")
species.coef$parameter <- trimws(species.coef$parameter)
species.coef$names <- trimws(species.coef$names)
head(species.coef)
species.coef <- species.coef[,c(1,5,6)]
str(species.coef)
species.coef <- species.coef %>%
pivot_wider(names_from = parameter, values_from = intrcpt)
species.coef <- as.data.frame(species.coef)
#correlation between lambda and lFC across species based on random effects from mod2
cor.test(species.coef$lambda, species.coef$lFC)
plot(species.coef$lambda, species.coef$lFC, main = "Model 2")
[[alternative HTML version deleted]]