Dear list, I'm quite new to MCMCglmm and after reading much of Jarrod Hadfield's documentation I was able to understand the basic construction and interpretation of models. However, I am having a bit of trouble when predicting values from my gaussian MCMC model. I am running bayesian phylogenetic mixed models to determine factors responsible for the variation of my movement parameter. My response variable is continuous and each point has a standard error associated to it and which I added as a statistical weight in the model using the mev argument in the MCMCglmm function. I have 5 random effects: phylogeny (to which I attributed the species covariance matrix using the ginverse argument), species, population, year and individual. My model is thus as follows: >>> model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration ???????????????????????????????????? , random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id ???????????????????????????????????? , family = "gaussian" ???????????????????????????????????? , mev = SE^2 # error variance associated to each data point ???????????????????????????????????? , ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo ???????????????????????????????????? , prior = prior1 ???????????????????????????????????? , data = Data ???????????????????????????????????? , nitt = 22e+04 ???????????????????????????????????? , burnin = 20000 ???????????????????????????????????? , thin = 100 ???????????????????????????????????? , pr=TRUE) When I try to predict values from the model above, I obtain the following error: >>> pred_DGLOB <- predict.MCMCglmm(model_D_GLOB_ch1 ?????????????????????????????????? , type = "response" ?????????????????????????????????? , interval = "confidence" ?????????????????????????????????? , newdata = newdt) >>> Error in rep(as.numeric(rcomponents %in% mcomponents), object$Random$nrt) : ????? invalid 'times' argument It appears that this error is originating from the following specific code in the predict.MCMCglmm function: >>> marginalise<-rep(as.numeric(rcomponents%in%mcomponents), object$Random$nrt) It appears that object$Random$nrt has one additional element compared to as.numeric(rcomponents%in%mcomponents). I already ran predict.MCMCglmm with models without any statistical weight and the predict.MCMCglmm function worked fine. Is this a problem caused by the fact that I have the standard error in the mev argument? If yes, is there a way to go around it? I may be skipping some steps before running the predict function but I did not find any documentation that addresses this issue. I would greatly appreciate your help and happy to provide further information if needed! Thank you in advance! Kamal
Predicting values from MCMCglmm model with statistical weight in mev argument
5 messages · Kamal Atmeh, Jarrod Hadfield
2 days later
Hi Kamal,
Can you post your sessionInfo()?
As a work around, use this model
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration
, random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units,
, family = "gaussian"
, ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo
, prior = prior1
, data = Data
, nitt = 22e+04
, burnin = 20000
, thin = 100
, pr=TRUE)
BUT make sure to fix the prior variance associated with the final random effect term (idh(SE):units) to one. Its identical to the model you've fitted, but the predict function should work.
Cheers,
Jarrod
On 15/02/2020 22:57, Kamal Atmeh wrote:
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration
, random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id
, family = "gaussian"
, mev = SE^2 # error variance associated to each data point
, ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo
, prior = prior1
, data = Data
, nitt = 22e+04
, burnin = 20000
, thin = 100
, pr=TRUE)
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
Hi Jarrod,
Thank you for your answer, the predict function worked! I used the
following non-informative prior with a fixed variance for the final
random effect as you suggested.
>>> prior1<-list(G=list(G1=list(V=1,nu=0.02)
??????????????????????? ,G2=list(V=1,nu=0.02)
??????????????????????? ,G3=list(V=1,nu=0.02)
??????????????????????? ,G4=list(V=1,nu=0.02)
??????????????????????? ,G5=list(V=1,nu=0.02)
??????????????????????? ,G6=list(V=1,fix=1)),
???????????????? R=list(V=1,nu=0.02))
model <- MCMCglmm(lD ~
tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration
???????????????????????????????????? , random =
~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units,
???????????????????????????????????? , family = "gaussian"
???????????????????????????????????? , ginverse = list(sp_phylo =
inv.phylo$Ainv) # include a custom matrix for argument phylo
???????????????????????????????????? , prior = prior1
???????????????????????????????????? , data = Data
???????????????????????????????????? , nitt = 22e+04
???????????????????????????????????? , burnin = 20000
???????????????????????????????????? , thin = 100
???????????????????????????????????? , pr=TRUE)
When doing expand.grid() to add in the predict function, I fixed the SE
parameter to the mean of all standard errors of my original data. Is
this a correct way to define the standard error column in my expand.grid
or should I choose one value as I did for the other random effects?
>>>>> newdt=expand.grid(tactic=c("F","H")
????????????????????? , period=c("PB","B")
????????????????????? , lbody=c(mean(Data$lbody),mean(Data$lbody) +
sd(Data$lbody))
????????????????????? , complique_KF=c("OU/OUF", "BM")
????????????????????? , mean.dhi_ndviqa_f.3=seq(min(Data
$mean.dhi_ndviqa_f.3), max(Data$mean.dhi_ndviqa_f.3),length.out = 500)
## When only hider, use length.out=500
????????????????????? , lintdur=c(mean(Data$lintdur),mean(Data$lintdur)
+ sd(Data$lintdur))
????????????????????? ,
lduration=c(mean(Data$lduration),mean(Data$lduration) + sd(Data$lduration))
????????????????????? , lnb.loc=c(mean(Data$lnb.loc),mean(Data
$lnb.loc) + sd(Data $lnb.loc))
????????????????????? , sp_phylo_glenn="Odo_hem"
????????????????????? , species2="Odo_hem"
????????????????????? , phylo_pop="Odo_hem-wyoming"
????????????????????? , phylo_popY="Ant_ame-red_desert-2015"
????????????????????? , phylo_pop_id="Bis_bis-PANP-1001"
? ? ?? >>>> ???? , SE=mean(Data$SE))? ## MEAN OF ALL STANDARD ERRORS IN
ORIGINAL DATA
I am posting below my sessionInfo() as you requested. Thanks again for
the help.
Cheers,
Kamal
R version 3.6.2 (2019-12-12)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)
Matrix products: default
locale:
[1] LC_COLLATE=French_France.1252? LC_CTYPE=French_France.1252
LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C?????????????????? LC_TIME=French_France.1252
attached base packages:
[1] grid????? parallel? stats???? graphics? grDevices utils datasets?
methods?? base
other attached packages:
?[1] PerformanceAnalytics_1.5.3 xts_0.11-2 zoo_1.8-6
?[4] plyr_1.8.4???????????????? SDMTools_1.1-221.1 ggthemes_4.2.0
?[7] SyncMove_0.1-0???????????? timeline_0.9 gtable_0.3.0
[10] plot3D_1.1.1?????????????? beepr_1.3 gatepoints_0.1.3
[13] RPostgreSQL_0.6-2????????? DBI_1.0.0 trajr_1.3.0
[16] scales_1.0.0?????????????? FactoMineR_1.42 factoextra_1.0.5
[19] dismo_1.1-4??????????????? raster_2.9-5 rgdal_1.4-4
[22] plotrix_3.7-6????????????? corrplot_0.84 adehabitatHR_0.4.16
[25] adehabitatLT_0.3.24??????? CircStats_0.2-6 boot_1.3-23
[28] adehabitatMA_0.3.13??????? deldir_0.1-22 maptools_0.9-5
[31] ks_1.11.5????????????????? influence.ME_0.9-9 visreg_2.5-1
[34] rgeos_0.4-3??????????????? sp_1.3-1 cowplot_0.9.4
[37] RColorBrewer_1.1-2???????? rgl_0.100.30 misc3d_0.8-4
[40] MCMCglmm_2.29????????????? coda_0.19-3 MASS_7.3-51.4
[43] adephylo_1.1-11??????????? egg_0.4.5 gridExtra_2.3
[46] plotly_4.9.0?????????????? ggplot2_3.2.0 phytools_0.6-99
[49] maps_3.3.0???????????????? ape_5.3 rptR_0.9.22
[52] sjPlot_2.8.2?????????????? nlme_3.1-142 ade4_1.7-13
[55] MuMIn_1.43.6?????????????? glmm_1.3.0 doParallel_1.0.14
[58] iterators_1.0.10?????????? foreach_1.4.4 mvtnorm_1.0-11
[61] trust_0.1-7??????????????? phylobase_0.8.6 lmerTest_3.1-0
[64] lme4_1.1-21??????????????? Matrix_1.2-18
loaded via a namespace (and not attached):
? [1] R.utils_2.9.0?????????? tidyselect_0.2.5 htmlwidgets_1.3
? [4] combinat_0.0-8????????? RNeXML_2.3.0 munsell_0.5.0
? [7] animation_2.6?????????? codetools_0.2-16 effectsize_0.1.1
?[10] units_0.6-3???????????? miniUI_0.1.1.1 withr_2.1.2
?[13] audio_0.1-6???????????? colorspace_1.4-1 knitr_1.23
?[16] uuid_0.1-2????????????? rstudioapi_0.10 leaps_3.0
?[19] stats4_3.6.2??????????? emmeans_1.4.4 mnormt_1.5-5
?[22] LearnBayes_2.15.1?????? vctrs_0.2.0 generics_0.0.2
?[25] clusterGeneration_1.3.4 xfun_0.8 itertools_0.1-3
?[28] adegenet_2.1.1????????? R6_2.4.0 manipulateWidget_0.10.0
?[31] assertthat_0.2.1??????? promises_1.0.1 phangorn_2.5.5
?[34] rlang_0.4.0???????????? zeallot_0.1.0 scatterplot3d_0.3-41
?[37] splines_3.6.2?????????? lazyeval_0.2.2 broom_0.5.2
?[40] reshape2_1.4.3????????? modelr_0.1.5 crosstalk_1.0.0
?[43] backports_1.1.4???????? httpuv_1.5.1 tensorA_0.36.1
?[46] tools_3.6.2???????????? spData_0.3.2 cubature_2.0.3
?[49] Rcpp_1.0.1????????????? progress_1.2.2 classInt_0.3-3
?[52] purrr_0.3.2???????????? prettyunits_1.0.2 haven_2.1.1
?[55] ggrepel_0.8.1?????????? cluster_2.1.0 magrittr_1.5
?[58] data.table_1.12.2?????? gmodels_2.18.1 sjmisc_2.8.3
?[61] hms_0.5.0?????????????? mime_0.7 xtable_1.8-4
?[64] XML_3.98-1.20?????????? sjstats_0.17.9 mclust_5.4.4
?[67] ggeffects_0.14.1??????? compiler_3.6.2 tibble_2.1.3
?[70] KernSmooth_2.23-16????? crayon_1.3.4 R.oo_1.22.0
?[73] minqa_1.2.4???????????? htmltools_0.3.6 mgcv_1.8-31
?[76] corpcor_1.6.9?????????? later_0.8.0 spdep_1.1-3
?[79] tidyr_1.0.2???????????? expm_0.999-4 sjlabelled_1.1.3
?[82] sf_0.7-6??????????????? permute_0.9-5 R.methodsS3_1.7.1
?[85] quadprog_1.5-7????????? gdata_2.18.0 insight_0.8.1
?[88] igraph_1.2.4.1????????? forcats_0.4.0 pkgconfig_2.0.2
?[91] flashClust_1.01-2?????? rncl_0.8.3 numDeriv_2016.8-1.1
?[94] foreign_0.8-72????????? xml2_1.2.1 webshot_0.5.1
?[97] estimability_1.3??????? stringr_1.4.0 digest_0.6.20
[100] parameters_0.5.0??????? vegan_2.5-6 fastmatch_1.1-0
[103] shiny_1.3.2???????????? gtools_3.8.1 nloptr_1.2.1
[106] lifecycle_0.1.0???????? jsonlite_1.6 seqinr_3.6-1
[109] viridisLite_0.3.0?????? pillar_1.4.2 lattice_0.20-38
[112] httr_1.4.0????????????? glue_1.3.1 bayestestR_0.5.2
[115] class_7.3-15??????????? stringi_1.4.3 performance_0.4.4
[118] dplyr_0.8.3???????????? e1071_1.7-2
Le 18/02/2020 ? 12:05, Jarrod Hadfield a ?crit?:
Hi Kamal, Can you post your sessionInfo()? As a work around, use this model model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration ???????????????????????????????????? , random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units, ???????????????????????????????????? , family = "gaussian" ???????????????????????????????????? , ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo ???????????????????????????????????? , prior = prior1 ???????????????????????????????????? , data = Data ???????????????????????????????????? , nitt = 22e+04 ???????????????????????????????????? , burnin = 20000 ???????????????????????????????????? , thin = 100 ???????????????????????????????????? , pr=TRUE) BUT make sure to fix the prior variance associated with the final random effect term (idh(SE):units) to one. Its identical to the model you've fitted, but the predict function should work. Cheers, Jarrod On 15/02/2020 22:57, Kamal Atmeh wrote:
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration ???????????????????????????????????? , random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id ???????????????????????????????????? , family = "gaussian" ???????????????????????????????????? , mev = SE^2 # error variance associated to each data point ???????????????????????????????????? , ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo ???????????????????????????????????? , prior = prior1 ???????????????????????????????????? , data = Data ???????????????????????????????????? , nitt = 22e+04 ???????????????????????????????????? , burnin = 20000 ???????????????????????????????????? , thin = 100 ???????????????????????????????????? , pr=TRUE)
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
Hi Kamal, It doesn?t matter what you set the SE to because you are a) marginalising the random effects b) the data are Gaussian and c) you are getting confidence intervals rather than prediction intervals. Cheers, Jarrod
On 18 Feb 2020, at 12:37, Kamal Atmeh <kamal.atmeh at hotmail.com<mailto:kamal.atmeh at hotmail.com>> wrote:
Hi Jarrod, Thank you for your answer, the predict function worked! I used the following non-informative prior with a fixed variance for the final random effect as you suggested.
prior1<-list(G=list(G1=list(V=1,nu=0.02)
,G2=list(V=1,nu=0.02)
,G3=list(V=1,nu=0.02)
,G4=list(V=1,nu=0.02)
,G5=list(V=1,nu=0.02)
,G6=list(V=1,fix=1)),
R=list(V=1,nu=0.02))
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration
, random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units,
, family = "gaussian"
, ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo
, prior = prior1
, data = Data
, nitt = 22e+04
, burnin = 20000
, thin = 100
, pr=TRUE)
When doing expand.grid() to add in the predict function, I fixed the SE parameter to the mean of all standard errors of my original data. Is this a correct way to define the standard error column in my expand.grid or should I choose one value as I did for the other random effects?
newdt=expand.grid(tactic=c("F","H")
, period=c("PB","B")
, lbody=c(mean(Data$lbody),mean(Data$lbody) + sd(Data$lbody))
, complique_KF=c("OU/OUF", "BM")
, mean.dhi_ndviqa_f.3=seq(min(Data $mean.dhi_ndviqa_f.3), max(Data$mean.dhi_ndviqa_f.3),length.out = 500) ## When only hider, use length.out=500
, lintdur=c(mean(Data$lintdur),mean(Data$lintdur) + sd(Data$lintdur))
, lduration=c(mean(Data$lduration),mean(Data$lduration) + sd(Data$lduration))
, lnb.loc=c(mean(Data$lnb.loc),mean(Data $lnb.loc) + sd(Data $lnb.loc))
, sp_phylo_glenn="Odo_hem"
, species2="Odo_hem"
, phylo_pop="Odo_hem-wyoming"
, phylo_popY="Ant_ame-red_desert-2015"
, phylo_pop_id="Bis_bis-PANP-1001"
>>>> , SE=mean(Data$SE)) ## MEAN OF ALL STANDARD ERRORS IN ORIGINAL DATA
I am posting below my sessionInfo() as you requested. Thanks again for the help.
Cheers,
Kamal
R version 3.6.2 (2019-12-12)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)
Matrix products: default
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252 LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C LC_TIME=French_France.1252
attached base packages:
[1] grid parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] PerformanceAnalytics_1.5.3 xts_0.11-2 zoo_1.8-6
[4] plyr_1.8.4 SDMTools_1.1-221.1 ggthemes_4.2.0
[7] SyncMove_0.1-0 timeline_0.9 gtable_0.3.0
[10] plot3D_1.1.1 beepr_1.3 gatepoints_0.1.3
[13] RPostgreSQL_0.6-2 DBI_1.0.0 trajr_1.3.0
[16] scales_1.0.0 FactoMineR_1.42 factoextra_1.0.5
[19] dismo_1.1-4 raster_2.9-5 rgdal_1.4-4
[22] plotrix_3.7-6 corrplot_0.84 adehabitatHR_0.4.16
[25] adehabitatLT_0.3.24 CircStats_0.2-6 boot_1.3-23
[28] adehabitatMA_0.3.13 deldir_0.1-22 maptools_0.9-5
[31] ks_1.11.5 influence.ME_0.9-9 visreg_2.5-1
[34] rgeos_0.4-3 sp_1.3-1 cowplot_0.9.4
[37] RColorBrewer_1.1-2 rgl_0.100.30 misc3d_0.8-4
[40] MCMCglmm_2.29 coda_0.19-3 MASS_7.3-51.4
[43] adephylo_1.1-11 egg_0.4.5 gridExtra_2.3
[46] plotly_4.9.0 ggplot2_3.2.0 phytools_0.6-99
[49] maps_3.3.0 ape_5.3 rptR_0.9.22
[52] sjPlot_2.8.2 nlme_3.1-142 ade4_1.7-13
[55] MuMIn_1.43.6 glmm_1.3.0 doParallel_1.0.14
[58] iterators_1.0.10 foreach_1.4.4 mvtnorm_1.0-11
[61] trust_0.1-7 phylobase_0.8.6 lmerTest_3.1-0
[64] lme4_1.1-21 Matrix_1.2-18
loaded via a namespace (and not attached):
[1] R.utils_2.9.0 tidyselect_0.2.5 htmlwidgets_1.3
[4] combinat_0.0-8 RNeXML_2.3.0 munsell_0.5.0
[7] animation_2.6 codetools_0.2-16 effectsize_0.1.1
[10] units_0.6-3 miniUI_0.1.1.1 withr_2.1.2
[13] audio_0.1-6 colorspace_1.4-1 knitr_1.23
[16] uuid_0.1-2 rstudioapi_0.10 leaps_3.0
[19] stats4_3.6.2 emmeans_1.4.4 mnormt_1.5-5
[22] LearnBayes_2.15.1 vctrs_0.2.0 generics_0.0.2
[25] clusterGeneration_1.3.4 xfun_0.8 itertools_0.1-3
[28] adegenet_2.1.1 R6_2.4.0 manipulateWidget_0.10.0
[31] assertthat_0.2.1 promises_1.0.1 phangorn_2.5.5
[34] rlang_0.4.0 zeallot_0.1.0 scatterplot3d_0.3-41
[37] splines_3.6.2 lazyeval_0.2.2 broom_0.5.2
[40] reshape2_1.4.3 modelr_0.1.5 crosstalk_1.0.0
[43] backports_1.1.4 httpuv_1.5.1 tensorA_0.36.1
[46] tools_3.6.2 spData_0.3.2 cubature_2.0.3
[49] Rcpp_1.0.1 progress_1.2.2 classInt_0.3-3
[52] purrr_0.3.2 prettyunits_1.0.2 haven_2.1.1
[55] ggrepel_0.8.1 cluster_2.1.0 magrittr_1.5
[58] data.table_1.12.2 gmodels_2.18.1 sjmisc_2.8.3
[61] hms_0.5.0 mime_0.7 xtable_1.8-4
[64] XML_3.98-1.20 sjstats_0.17.9 mclust_5.4.4
[67] ggeffects_0.14.1 compiler_3.6.2 tibble_2.1.3
[70] KernSmooth_2.23-16 crayon_1.3.4 R.oo_1.22.0
[73] minqa_1.2.4 htmltools_0.3.6 mgcv_1.8-31
[76] corpcor_1.6.9 later_0.8.0 spdep_1.1-3
[79] tidyr_1.0.2 expm_0.999-4 sjlabelled_1.1.3
[82] sf_0.7-6 permute_0.9-5 R.methodsS3_1.7.1
[85] quadprog_1.5-7 gdata_2.18.0 insight_0.8.1
[88] igraph_1.2.4.1 forcats_0.4.0 pkgconfig_2.0.2
[91] flashClust_1.01-2 rncl_0.8.3 numDeriv_2016.8-1.1
[94] foreign_0.8-72 xml2_1.2.1 webshot_0.5.1
[97] estimability_1.3 stringr_1.4.0 digest_0.6.20
[100] parameters_0.5.0 vegan_2.5-6 fastmatch_1.1-0
[103] shiny_1.3.2 gtools_3.8.1 nloptr_1.2.1
[106] lifecycle_0.1.0 jsonlite_1.6 seqinr_3.6-1
[109] viridisLite_0.3.0 pillar_1.4.2 lattice_0.20-38
[112] httr_1.4.0 glue_1.3.1 bayestestR_0.5.2
[115] class_7.3-15 stringi_1.4.3 performance_0.4.4
[118] dplyr_0.8.3 e1071_1.7-2
Le 18/02/2020 ? 12:05, Jarrod Hadfield a ?crit :
Hi Kamal,
Can you post your sessionInfo()?
As a work around, use this model
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration
, random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units,
, family = "gaussian"
, ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo
, prior = prior1
, data = Data
, nitt = 22e+04
, burnin = 20000
, thin = 100
, pr=TRUE)
BUT make sure to fix the prior variance associated with the final random effect term (idh(SE):units) to one. Its identical to the model you've fitted, but the predict function should work.
Cheers,
Jarrod
On 15/02/2020 22:57, Kamal Atmeh wrote:
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration
, random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id
, family = "gaussian"
, mev = SE^2 # error variance associated to each data point
, ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo
, prior = prior1
, data = Data
, nitt = 22e+04
, burnin = 20000
, thin = 100
, pr=TRUE)
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
Hi Jarrod, Perfect. Thank you very much for your answer and for your help. Cheers, Kamal Le 18/02/2020 ? 19:25, HADFIELD Jarrod a ?crit?:
Hi Kamal, It doesn?t matter what you set the SE to because you are a) marginalising the random effects b) the data are Gaussian and c) you are getting confidence intervals rather than prediction intervals. Cheers, Jarrod
On 18 Feb 2020, at 12:37, Kamal Atmeh <kamal.atmeh at hotmail.com <mailto:kamal.atmeh at hotmail.com>> wrote: Hi Jarrod, Thank you for your answer, the predict function worked! I used the following non-informative prior with a fixed variance for the final random effect as you suggested.
prior1<-list(G=list(G1=list(V=1,nu=0.02)
??????????????????????? ,G2=list(V=1,nu=0.02) ??????????????????????? ,G3=list(V=1,nu=0.02) ??????????????????????? ,G4=list(V=1,nu=0.02) ??????????????????????? ,G5=list(V=1,nu=0.02) ??????????????????????? ,G6=list(V=1,fix=1)), ???????????????? R=list(V=1,nu=0.02)) model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration ???????????????????????????????????? , random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units, ???????????????????????????????????? , family = "gaussian" ???????????????????????????????????? , ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo ???????????????????????????????????? , prior = prior1 ???????????????????????????????????? , data = Data ???????????????????????????????????? , nitt = 22e+04 ???????????????????????????????????? , burnin = 20000 ???????????????????????????????????? , thin = 100 ???????????????????????????????????? , pr=TRUE) When doing expand.grid() to add in the predict function, I fixed the SE parameter to the mean of all standard errors of my original data. Is this a correct way to define the standard error column in my expand.grid or should I choose one value as I did for the other random effects?
newdt=expand.grid(tactic=c("F","H")
????????????????????? , period=c("PB","B")
????????????????????? , lbody=c(mean(Data$lbody),mean(Data$lbody) +
sd(Data$lbody))
????????????????????? , complique_KF=c("OU/OUF", "BM")
????????????????????? , mean.dhi_ndviqa_f.3=seq(min(Data
$mean.dhi_ndviqa_f.3), max(Data$mean.dhi_ndviqa_f.3),length.out =
500) ## When only hider, use length.out=500
????????????????????? ,
lintdur=c(mean(Data$lintdur),mean(Data$lintdur) + sd(Data$lintdur))
????????????????????? ,
lduration=c(mean(Data$lduration),mean(Data$lduration) +
sd(Data$lduration))
????????????????????? , lnb.loc=c(mean(Data$lnb.loc),mean(Data
$lnb.loc) + sd(Data $lnb.loc))
????????????????????? , sp_phylo_glenn="Odo_hem"
????????????????????? , species2="Odo_hem"
????????????????????? , phylo_pop="Odo_hem-wyoming"
????????????????????? , phylo_popY="Ant_ame-red_desert-2015"
????????????????????? , phylo_pop_id="Bis_bis-PANP-1001"
? ? ?? >>>> ???? , SE=mean(Data$SE))? ## MEAN OF ALL STANDARD ERRORS
IN ORIGINAL DATA
I am posting below my sessionInfo() as you requested. Thanks again
for the help.
Cheers,
Kamal
R version 3.6.2 (2019-12-12)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)
Matrix products: default
locale:
[1] LC_COLLATE=French_France.1252 LC_CTYPE=French_France.1252
LC_MONETARY=French_France.1252
[4] LC_NUMERIC=C LC_TIME=French_France.1252
attached base packages:
[1] grid????? parallel? stats???? graphics? grDevices utils????
datasets? methods?? base
other attached packages:
?[1] PerformanceAnalytics_1.5.3 xts_0.11-2???????????????? zoo_1.8-6
?[4] plyr_1.8.4 SDMTools_1.1-221.1???????? ggthemes_4.2.0
?[7] SyncMove_0.1-0 timeline_0.9?????????????? gtable_0.3.0
[10] plot3D_1.1.1 beepr_1.3????????????????? gatepoints_0.1.3
[13] RPostgreSQL_0.6-2 DBI_1.0.0????????????????? trajr_1.3.0
[16] scales_1.0.0 FactoMineR_1.42??????????? factoextra_1.0.5
[19] dismo_1.1-4 raster_2.9-5?????????????? rgdal_1.4-4
[22] plotrix_3.7-6 corrplot_0.84????????????? adehabitatHR_0.4.16
[25] adehabitatLT_0.3.24 CircStats_0.2-6??????????? boot_1.3-23
[28] adehabitatMA_0.3.13 deldir_0.1-22????????????? maptools_0.9-5
[31] ks_1.11.5 influence.ME_0.9-9???????? visreg_2.5-1
[34] rgeos_0.4-3 sp_1.3-1?????????????????? cowplot_0.9.4
[37] RColorBrewer_1.1-2 rgl_0.100.30?????????????? misc3d_0.8-4
[40] MCMCglmm_2.29 coda_0.19-3??????????????? MASS_7.3-51.4
[43] adephylo_1.1-11 egg_0.4.5????????????????? gridExtra_2.3
[46] plotly_4.9.0 ggplot2_3.2.0????????????? phytools_0.6-99
[49] maps_3.3.0 ape_5.3??????????????????? rptR_0.9.22
[52] sjPlot_2.8.2 nlme_3.1-142?????????????? ade4_1.7-13
[55] MuMIn_1.43.6 glmm_1.3.0???????????????? doParallel_1.0.14
[58] iterators_1.0.10 foreach_1.4.4????????????? mvtnorm_1.0-11
[61] trust_0.1-7 phylobase_0.8.6??????????? lmerTest_3.1-0
[64] lme4_1.1-21 Matrix_1.2-18
loaded via a namespace (and not attached):
? [1] R.utils_2.9.0?????????? tidyselect_0.2.5 htmlwidgets_1.3
? [4] combinat_0.0-8????????? RNeXML_2.3.0 munsell_0.5.0
? [7] animation_2.6?????????? codetools_0.2-16 effectsize_0.1.1
?[10] units_0.6-3???????????? miniUI_0.1.1.1 withr_2.1.2
?[13] audio_0.1-6???????????? colorspace_1.4-1 knitr_1.23
?[16] uuid_0.1-2????????????? rstudioapi_0.10 leaps_3.0
?[19] stats4_3.6.2??????????? emmeans_1.4.4 mnormt_1.5-5
?[22] LearnBayes_2.15.1?????? vctrs_0.2.0 generics_0.0.2
?[25] clusterGeneration_1.3.4 xfun_0.8 itertools_0.1-3
?[28] adegenet_2.1.1????????? R6_2.4.0 manipulateWidget_0.10.0
?[31] assertthat_0.2.1??????? promises_1.0.1 phangorn_2.5.5
?[34] rlang_0.4.0???????????? zeallot_0.1.0 scatterplot3d_0.3-41
?[37] splines_3.6.2?????????? lazyeval_0.2.2 broom_0.5.2
?[40] reshape2_1.4.3????????? modelr_0.1.5 crosstalk_1.0.0
?[43] backports_1.1.4???????? httpuv_1.5.1 tensorA_0.36.1
?[46] tools_3.6.2???????????? spData_0.3.2 cubature_2.0.3
?[49] Rcpp_1.0.1????????????? progress_1.2.2 classInt_0.3-3
?[52] purrr_0.3.2???????????? prettyunits_1.0.2 haven_2.1.1
?[55] ggrepel_0.8.1?????????? cluster_2.1.0 magrittr_1.5
?[58] data.table_1.12.2?????? gmodels_2.18.1 sjmisc_2.8.3
?[61] hms_0.5.0?????????????? mime_0.7 xtable_1.8-4
?[64] XML_3.98-1.20?????????? sjstats_0.17.9 mclust_5.4.4
?[67] ggeffects_0.14.1??????? compiler_3.6.2 tibble_2.1.3
?[70] KernSmooth_2.23-16????? crayon_1.3.4 R.oo_1.22.0
?[73] minqa_1.2.4???????????? htmltools_0.3.6 mgcv_1.8-31
?[76] corpcor_1.6.9?????????? later_0.8.0 spdep_1.1-3
?[79] tidyr_1.0.2???????????? expm_0.999-4 sjlabelled_1.1.3
?[82] sf_0.7-6??????????????? permute_0.9-5 R.methodsS3_1.7.1
?[85] quadprog_1.5-7????????? gdata_2.18.0 insight_0.8.1
?[88] igraph_1.2.4.1????????? forcats_0.4.0 pkgconfig_2.0.2
?[91] flashClust_1.01-2?????? rncl_0.8.3 numDeriv_2016.8-1.1
?[94] foreign_0.8-72????????? xml2_1.2.1 webshot_0.5.1
?[97] estimability_1.3??????? stringr_1.4.0 digest_0.6.20
[100] parameters_0.5.0??????? vegan_2.5-6 fastmatch_1.1-0
[103] shiny_1.3.2???????????? gtools_3.8.1 nloptr_1.2.1
[106] lifecycle_0.1.0???????? jsonlite_1.6 seqinr_3.6-1
[109] viridisLite_0.3.0?????? pillar_1.4.2 lattice_0.20-38
[112] httr_1.4.0????????????? glue_1.3.1 bayestestR_0.5.2
[115] class_7.3-15??????????? stringi_1.4.3 performance_0.4.4
[118] dplyr_0.8.3???????????? e1071_1.7-2
Le 18/02/2020 ? 12:05, Jarrod Hadfield a ?crit?:
Hi Kamal, Can you post your sessionInfo()? As a work around, use this model model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration ???????????????????????????????????? , random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id + idh(SE):units, ???????????????????????????????????? , family = "gaussian" ???????????????????????????????????? , ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo ???????????????????????????????????? , prior = prior1 ???????????????????????????????????? , data = Data ???????????????????????????????????? , nitt = 22e+04 ???????????????????????????????????? , burnin = 20000 ???????????????????????????????????? , thin = 100 ???????????????????????????????????? , pr=TRUE) BUT make sure to fix the prior variance associated with the final random effect term (idh(SE):units) to one. Its identical to the model you've fitted, but the predict function should work. Cheers, Jarrod On 15/02/2020 22:57, Kamal Atmeh wrote:
model <- MCMCglmm(lD ~ tactic*period*seasonality+complique_KF+lbody+lintdur+lnb.loc+lduration ???????????????????????????????????? , random = ~sp_phylo+species2+phylo_pop+phylo_popY+phylo_pop_id ???????????????????????????????????? , family = "gaussian" ???????????????????????????????????? , mev = SE^2 # error variance associated to each data point ???????????????????????????????????? , ginverse = list(sp_phylo = inv.phylo$Ainv) # include a custom matrix for argument phylo ???????????????????????????????????? , prior = prior1 ???????????????????????????????????? , data = Data ???????????????????????????????????? , nitt = 22e+04 ???????????????????????????????????? , burnin = 20000 ???????????????????????????????????? , thin = 100 ???????????????????????????????????? , pr=TRUE)
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.