search

willpearse / pez

Phylogenetics for the Environmental Sciences
GNU General Public License v3.0
14 stars 7 forks source link

Error in random effects-"Number of items to replace is not a multiple of replacement length" #72

Open Jigyasa3 opened 4 years ago

Jigyasa3 commented 4 years ago

Hey!

I am following the supplementary material of [https://nph.onlinelibrary.wiley.com/doi/epdf/10.1111/nph.14397] on my own data. My data is slightly different, but the questions I am looking at are similar to the paper. But I keep running into an error, which I can't find any information about on the internet.

My aim is to see if the lineage of the species helps explain the gene relative abundance.

my data-

K00297 is the log transformed gene relative abundance data-

> wood_0.5_df_sub_pez2%>%dplyr::select(K00297,lineage,diet_type_number2)%>%head() K00297 lineage diet_type_number2 tip_names 1 0.4009284 Termitinae 3 Termitinae_sp1 2 0.3649912 Termitinae 3 Termitinae_sp2 3 -0.1522581 Termitinae 3 Termitinae_sp3 4 0.1961155 Termitinae 2 Termitinae_sp4 5 0.5761659 Termitinae 2 Termitinae_sp5 6 0.2031437 Termitinae 3 Termitinae_sp6

the R code-

>wood_0.5_df_sub_pez2$lineage = as.factor(wood_0.5_df_sub_pez2$lineage) 
>wood_0.5_df_sub_pez2$diet_type_number2 = as.factor(wood_0.5_df_sub_pez2$diet_type_number2) 
>nspp <- nlevels(wood_0.5_df_sub_pez2$lineage)
 >nsite <- nlevels(wood_0.5_df_sub_pez2$diet_type_number2)

to standardize the phylo var-covar matrix

>Vphy <- vcv(phy)
>Vphy <- Vphy[order(phy$tip.label),order(phy$tip.label)]
>Vphy <- Vphy / max(Vphy)
>Vphy <- Vphy / det(Vphy) ^ (1 / nspp)

specify random effects

re.site <- list(1, site = wood_0.5_df_sub_pez2$diet_type_number2, covar = diag(nsite))
  re.sp <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = diag(nspp))
  re.sp.phy <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = Vphy)

sp is nested within site

re.nested.phy <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = Vphy, site = wood_0.5_df_sub_pez2$diet_type_number2)
 re.nested.rep <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = solve(Vphy), site = wood_0.5_df_sub_pez2$diet_type_number2)

run the full model without s2.init variable to get the random effect variance, and then add that to the final model.

this is the full model,without s2.init

z_test <- communityPGLMM(K00297 ~ 1, data = wood_0.5_df_sub_pez2, family = "gaussian", sp = wood_0.5_df_sub_pez2$lineage, site=wood_0.5_df_sub_pez2$diet_type_number2, random.effects = list(re.sp, re.sp.phy,re.site, re.nested.phy), REML = F, verbose = F)

the error-

Error in Z.i[, counter] <- re.i[[1]] * as.numeric(i.levels == re.i[[2]]) : number of items to replace is not a multiple of replacement length

If I run the model with sp=wood_0.5_df_sub_pez2$tip.names (and make the corresponding changes to the random effects also), then the error is-

Error: Matrices must have same dimensions in .Arith.Csparse(e1, e2, .Generic, class. = "dgCMatrix")

Any suggestions? Thanks!

willpearse commented 4 years ago

Hello,

Thanks for this. It sounds like you've got a problem whereby the dimensions of your data don't match, and so my advice would be to double-check those dimensions and, if you can't see anything obvious there, use the comparative.comm data format within pez to structure your data, then use as.data.frame to turn that structure into something commnuityPGLMM can work with. You'll then find this error goes away, I think.

Let me know if that helps. If it doesn't, my post-doc and I have made a series of tutorials for how to run PGLMM using rstan; if that's something that's helpful to you, I can send those to you if you email me at will.pearse@usu.edu (I am about to go away travelling for a few days though, so I may be slow to reply). I would encourage you to try the as.data.frame method first, though.

Cheers,

Will

Measuring phylogenetic structure? Try install.packages('pez') I have a baby and so may be slow to reply to emails. Thank you for your patience.

Will Pearse http://www.pearselab.com/ Assistant Professor of Biology, Utah State University Office: +1-435-797-0831; Room LSB-333 Skype: will.pearse

On Fri, 10 Jan 2020 at 01:09, Jigyasa3 notifications@github.com wrote:

Hey!

I am following the supplementary material of [ https://nph.onlinelibrary.wiley.com/doi/epdf/10.1111/nph.14397] on my own data. My data is slightly different, but the questions I am looking at are similar to the paper. But I keep running into an error, which I can't find any information about on the internet.

My aim is to see if the lineage of the species helps explain the gene relative abundance.

my data-

K00297 is the log transformed gene relative abundance data-

wood_0.5_df_sub_pez2%>%dplyr::select(K00297,lineage,diet_type_number2)%>%head() K00297 lineage diet_type_number2 tip_names 1 0.4009284 Termitinae 3 Termitinae_sp1 2 0.3649912 Termitinae 3 Termitinae_sp2 3 -0.1522581 Termitinae 3 Termitinae_sp3 4 0.1961155 Termitinae 2 Termitinae_sp4 5 0.5761659 Termitinae 2 Termitinae_sp5 6 0.2031437 Termitinae 3 Termitinae_sp6

the R code-

wood_0.5_df_sub_pez2$lineage = as.factor(wood_0.5_df_sub_pez2$lineage) wood_0.5_df_sub_pez2$diet_type_number2 = as.factor(wood_0.5_df_sub_pez2$diet_type_number2) nspp <- nlevels(wood_0.5_df_sub_pez2$lineage) nsite <- nlevels(wood_0.5_df_sub_pez2$diet_type_number2)

to standardize the phylo var-covar matrix

Vphy <- vcv(phy) Vphy <- Vphy[order(phy$tip.label),order(phy$tip.label)] Vphy <- Vphy / max(Vphy) Vphy <- Vphy / det(Vphy) ^ (1 / nspp)

specify random effects

re.site <- list(1, site = wood_0.5_df_sub_pez2$diet_type_number2, covar = diag(nsite)) re.sp <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = diag(nspp)) re.sp.phy <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = Vphy)

sp is nested within site

re.nested.phy <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = Vphy, site = wood_0.5_df_sub_pez2$diet_type_number2) re.nested.rep <- list(1, sp = wood_0.5_df_sub_pez2$lineage, covar = solve(Vphy), site = wood_0.5_df_sub_pez2$diet_type_number2)

run the full model without s2.init variable to get the random effect

variance, and then add that to the final model. this is the full model,without s2.init

z_test <- communityPGLMM(K00297 ~ 1, data = wood_0.5_df_sub_pez2, family = "gaussian", sp = wood_0.5_df_sub_pez2$lineage, site=wood_0.5_df_sub_pez2$diet_type_number2, random.effects = list(re.sp, re.sp.phy,re.site, re.nested.phy), REML = F, verbose = F)

the error-

Error in Z.i[, counter] <- re.i[[1]] * as.numeric(i.levels == re.i[[2]]) : number of items to replace is not a multiple of replacement length

If I run the model with sp=wood_0.5_df_sub_pez2$tip.names (and make the

corresponding changes to the random effects also), then the error is- Error: Matrices must have same dimensions in .Arith.Csparse(e1, e2, .Generic, class. = "dgCMatrix")

Any suggestions? Thanks!

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/willpearse/pez/issues/72?email_source=notifications&email_token=AAJNYUQEI2Y3BMUEKPEOTSLQ5AUKBA5CNFSM4KFD7MT2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4IFIXCAA, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJNYUT7VHYCTALCD26FEADQ5AUKBANCNFSM4KFD7MTQ .