[R-sig-ME] Predicting values from MCMCglmm model with statistical weight in mev argument
HADFIELD Jarrod
j@h@d||e|d @end|ng |rom ed@@c@uk
Tue Feb 18 19:25:58 CET 2020
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 using hotmail.com<mailto:kamal.atmeh using 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.
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list