[R-sig-ME] Predicting values from MCMCglmm model with statistical weight in mev argument
Kamal Atmeh
k@m@|@@tmeh @end|ng |rom hotm@||@com
Wed Feb 19 15:00:34 CET 2020
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 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