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

[Jigyasa3]

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]

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 ***@***.***> 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 <#72?email_source=notifications&email_token=AAJNYUQEI2Y3BMUEKPEOTSLQ5AUKBA5CNFSM4KFD7MT2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4IFIXCAA>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAJNYUT7VHYCTALCD26FEADQ5AUKBANCNFSM4KFD7MTQ> .