murphymv / semEff

Automatic Calculation of Effects for Piecewise Structural Equation Models
https://murphymv.github.io/semEff/
GNU General Public License v3.0
10 stars 0 forks source link

SemEff using gls models #53

Open carvazgon90 opened 9 months ago

carvazgon90 commented 9 months ago

Hello!

First of all, thank you for creating this wonderful R package 😊.

I have fitted the following psem, which consists of a list of phylogenetic least squared models (gls):

psem(gls(genome_size ~ temperature + precipitation, tree.y, correlation = corBrownian(phy= newT2)), gls(PC1_def ~ genome_size + temperature + precipitation, data.y, correlation = corBrownian(phy= newT2)), gls(PC2_def ~ genome_size + temperature + precipitation, data.y, correlation = corBrownian(phy= newT2)), gls(PC3_def ~ genome_size + temperature + precipitation, data.y, correlation = corBrownian(phy= newT2)), gls(herbivory ~ PC1_def + PC2_def + PC3_def + genome_size + temperature + precipitation, data.y, correlation = corBrownian(phy= newT2)))

However, when I tried to obtain total, direct, and indirect effects using 'semEff()', I encountered the following error:

Error in getData(m) : 'data' does not contain all variables used to fit model.

I assume this error occurs because 'gls' models require both the data and a correlation structure based on a phylogenetic tree, which is not included in the 'data' object.

Is there any way to resolve this issue? I would greatly appreciate your help.

Thanks!

Carla

murphymv commented 9 months ago

Hi Carla,

Yes it looks like you're right. For bootstrapping, semEff requires all variables to be in 'data'. I could have a look at the code to see if I can include the correlation structure in the data prior to resampling. I'll let you know.

Thanks for highlighting the issue.

Cheers, Mark

murphymv commented 9 months ago

Is the variable 'newT2' in 'data.y'? If so, I don't think there should be an issue as the function getData() does actually evaluate whatever's in the correlation argument of gls() (having looked at the code).

carvazgon90 commented 9 months ago

Dear Mark,

Thank you very much for your help.

When fitting PSEMs using the gls() function, the phylogenetic tree (newT2) is not within the same data object; instead, it exists as a separate object.

I also attempted to fit PSEMs using another function for phylogenetic models, the pgls() function, which takes a combined object as input, including both data and the phylogenetic tree. However, when I try to run semEff I get the following error:

Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) : cannot coerce class ‘"comparative.data"’ to a data.frame

The code for psems with pgls is as follows (using piecewiseSEM v. 1.2.1)

phylogeny.y <- comparative.data(newT2, data.y, SP, vcv=TRUE, vcv.dim=3)

list.pgls.y = list( pgls(genome_size ~ temperature + precipitation, phylogeny.y, lambda='ML'), pgls(PC1_def ~ genome_size + temperature + precipitation, phylogeny.y, lambda='ML'), pgls(PC2_def ~ genome_size + temperature + precipitation, phylogeny.y, lambda='ML'), pgls(PC3_def ~ genome_size + temperature + precipitation, phylogeny.y, lambda='ML'), pgls(herbivory ~ PC1_def + PC2_def + PC3_def + genome_size + temperature + precipitation, phylogeny.y, lambda='ML'))

psem.pgls.y = sem.fit(list.pgls.y, phylogeny.y)

Thanks again for you help,

respectfully,

Carla

Carla Vázquez-González Postdoctoral researcher (GAIN - Fulbright)

twitter.com/CarlaVzquezGon1http://twitter.com/CarlaVzquezGon1 researchgate.net/profile/Carla-Vazquez-Gonzalezhttp://researchgate.net/profile/Carla-Vazquez-Gonzalez

*Department of Ecology and Evolutionary Biology - Mooney's Lab University of California Irvine. CA. USA.

https://tritrophic.weebly.com/

*Evolutionary Ecology of Plant-Herbivore Interactions Misión Biológica de Galicia (CSIC) Galicia, Spain www.plantherbivory.weebly.comhttp://www.plantherbivory.weebly.com/


De: Mark Murphy @.> Enviado: domingo, 1 de octubre de 2023 10:39 Para: murphymv/semEff @.> Cc: Carla Vazquez Gonzalez @.>; Author @.> Asunto: Re: [murphymv/semEff] SemEff using gls models (Issue #53)

Is the variable 'newT2' in 'data.y'? If so, I don't think there should be an issue as the function getData() evaluates whatever's in the correlation argument of gls().

— Reply to this email directly, view it on GitHubhttps://urldefense.com/v3/__https://github.com/murphymv/semEff/issues/53*issuecomment-1742148957__;Iw!!CzAuKJ42GuquVTTmVmPViYEvSg!NmZPZ_4Hsue9Mp7buHUbcNYaGKmtva6Re_YD9yRN0eDNJAhh0su5K_rT98X8wrbr6Sg4yAJ_eSWcjFrVLA2-TaHF$, or unsubscribehttps://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/BC2ENRPZBMHOAHETJI4GP5LX5GTEFANCNFSM6AAAAAA5GVRJ7A__;!!CzAuKJ42GuquVTTmVmPViYEvSg!NmZPZ_4Hsue9Mp7buHUbcNYaGKmtva6Re_YD9yRN0eDNJAhh0su5K_rT98X8wrbr6Sg4yAJ_eSWcjFrVLFb6JgmX$. You are receiving this because you authored the thread.Message ID: @.***>

dumaskvn commented 3 months ago

Good evening dear all, I'm upping this issue as I face it as well, and I was wondering whether a solution had been found ?

So far I traced the Error in getData(m) : 'data' does not contain all variables used to fit model. error to line 105 of stdEff-fun.R file, were the "vn" variable contain the name of the phyogenetic tree used to specify the correlation, which is unlikely to be found in the model data.frame data. Removing the name of this tree from "vn" solved the issue and made the "getData" function to work.

my modifictation : line 105 : if (!all(!vn[grepl("tree.name",vn)] %in% names(d)))

It seems to me that the other functions of this package then rely on the "update" function to compute metrics on models, which function replicate the correlation structure specified at the origin in gls models. However, I can't get the bootstrapping function to work on gls models, and my understanding of the functions is limited..

line 352 of bootEff-fun.R BE <- rMapply(bootEff, m, w, SIMPLIFY = FALSE)

rMapply(bootEff, m, w, SIMPLIFY = FALSE)

ORDINARY NONPARAMETRIC BOOTSTRAP

Call: boot::boot(data = x, statistic = s, R = R, parallel = parallel, ncpus = nc, cl = cl)

Bootstrap Statistics : WARNING: All values of t1* are NA Warning messages: 1: In (function (m, w) : 101 model fit(s) or parameter estimation(s) failed. NAs reported/generated. 2: In attributes(B$t)[c("sim", "seed", "n")] <- c(B$sim, seed, nrow(B$data)) : number of items to replace is not a multiple of replacement length

Could you explain what is going on here ?

Sincerely yours, Keyvan