Message-ID: <DB7PR07MB507721EFE3925FF445CD0FBDFE110@DB7PR07MB5077.eurprd07.prod.outlook.com>
Date: 2020-02-18T12:37:29Z
From: Kamal Atmeh
Subject: Predicting values from MCMCglmm model with statistical weight in mev argument
In-Reply-To: <54f40b10-e31b-6310-e93e-f6ea1c4646da@ed.ac.uk>
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.
[[alternative HTML version deleted]]