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.