diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..d821302 --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,4 @@ +^renv$ +^renv\.lock$ +^.*\.Rproj$ +^\.Rproj\.user$ diff --git a/.Rprofile b/.Rprofile new file mode 100644 index 0000000..e35a681 --- /dev/null +++ b/.Rprofile @@ -0,0 +1,3 @@ +source("renv/activate.R") + +devtools::load_all() diff --git a/.drake/.gitignore b/.drake/.gitignore new file mode 100644 index 0000000..120f485 --- /dev/null +++ b/.drake/.gitignore @@ -0,0 +1,2 @@ +* +!/.gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fbbf5da --- /dev/null +++ b/.gitignore @@ -0,0 +1,21 @@ +# R auxiliary files +.Rproj.user +.Rhistory +.RData +.Ruserdata + +# Other files types + +# htlm output +Todo.html +Documents/*.html +*.docx + +# Project folders + +# data_raw +*.rda +*.rds + +# Mac auxiliary files +.DS_Store diff --git a/Analysis/Analysis_Drake.R b/Analysis/Analysis_Drake.R new file mode 100644 index 0000000..2e4ea67 --- /dev/null +++ b/Analysis/Analysis_Drake.R @@ -0,0 +1,84 @@ +######################## +#### Analysis #### +######################## + +#---- load env ---- +env <- devtools::load_all()$env + +#---- Plan ----- + +# load plan +plan <- get_analysis_plan() + +# Configure the analysis plan +config <- drake::drake_config(plan, + prework = "devtools::load_all()", + envir = env) + +# Plot the analysis plan +drake::vis_drake_graph(config, font_size = 16, targets_only = FALSE) + +#---- Make ---- + +# Delate the analysis results +# drake::clean(destroy = TRUE) + +# Run the analysis +drake::make(prework = "devtools::load_all()", + envir = env, + plan = plan) + +# Plot the analysis plan +drake::vis_drake_graph(config, font_size = 16, targets_only = FALSE) + +#---- load ---- + +# data +drake::loadd(data_munged) + +#---- cluster ---- + +# cluster tat +drake::loadd(cluster_mother_tat) +drake::loadd(cluster_father_tat) + +plot(cluster_mother_tat, main = "Dendrogramma") +rect.hclust(cluster_mother_tat, k=4, border="red") +plot(cluster_father_tat, main = "Dendrogramma") +rect.hclust(cluster_father_tat, k=4, border="red") + +drake::loadd(data_cluster_tat) +table(data_cluster_tat$cluster_mother_tat) +table(data_cluster_tat$cluster_father_tat) + +plot_scores_mother(data = data_cluster_tat) +plot_scores_father(data = data_cluster_tat) + +# mclust +drake::loadd(mclust_mother) +drake::loadd(mclust_father) + +plot(mclust_mother) +summary(mclust_mother) + +plot(mclust_father) +summary(mclust_father) + +#---- brms models ---- + +drake::loadd(fit_int) + +summary(fit_int) +plot(brms::conditional_effects(fit_int), ask = FALSE) + +plot(fit_int) + +drake::loadd(stan_data) + +str(fit_int$fit) +#----- + + + + +#---- diff --git a/Analysis/draft.R b/Analysis/draft.R new file mode 100644 index 0000000..67463f7 --- /dev/null +++ b/Analysis/draft.R @@ -0,0 +1,115 @@ +load("Data/ddef_datiECR.rda") +ddef=ddef_datiECR[ddef_datiECR$age_year<12.30,] +ddef<- subset(ddef, ddef$fas !="NA") +nrow(ddef) #847 + +#---- cluster madre ---- + +dcM <- dist(ddef[,c(60:71)], method = "euclidean") # distance matrix +fit <- hclust(dcM, method="ward.D") +plot(fit, main = "Dendrogramma") +# dal dendrogramma sembrano emergere 4 cluster +rect.hclust(fit, k=4, border="red") +ddef$gruppiattM <- as.factor(cutree(fit, k=4)) +table(ddef$gruppiattM) + + +# profilo +fattore=rep(c("Ansia","Evitamento"),rep(dim(ddef)[1],2)) +y=c(scale(ddef$Anxm),scale(ddef$Avm)) +y1=c(ddef$Anxm,ddef$Avm) +gruppo=rep(ddef$gruppiattM,2) +dprof=data.frame(y,y1,fattore,gruppo) + + +dprof$nomi.gruppi=factor(dprof$gruppo,labels=c("Gruppo1: Sicuri (n = 231)","Gruppo 2: Ansiosi (n = 286)","Gruppo 3: Ansiosi-Evitanti (n = 100)","Gruppo 4: Evitanti (n = 230)")) +par(mar=c(6, 5, 4, 2) + 0.1) + +sciplot::bargraph.CI(fattore, y, group = nomi.gruppi, data=dprof,legend=T,ylim=c(-2.4,1.6),x.leg=1.2, y.leg=-1.2,ylab="Punteggi standardizzati",xlab="",cex.lab=.8,cex.leg=.8, col=c("yellow","green","blue","red")) + + +#---- cluster padre ---- + +dcP <- dist(ddef[,c(99:110)], method = "euclidean") # distance matrix +fit <- hclust(dcP, method="ward.D") +plot(fit, main = "Dendrogramma") +# dal dendrogramma sembrano emergere 4 cluster +rect.hclust(fit, k=4, border="red") +ddef$gruppiattP <- as.factor(cutree(fit, k=4)) + +# Profilo +fattore=rep(c("Ansia","Evitamento"),rep(dim(ddef)[1],2)) +y=c(scale(ddef$Anxp),scale(ddef$Avp)) +y1=c(ddef$Anxp,ddef$Avp) +gruppo=rep(ddef$gruppiattP,2) +dprof=data.frame(y,y1,fattore,gruppo) + +dprof$nomi.gruppi=factor(dprof$gruppo,labels=c("Gruppo1: Sicuri (n = 206)","Gruppo 2: Ansiosi (n = 230)","Gruppo 3: Evitanti (n = 311 )","Gruppo 4: Ansiosi-Evitanti (n = 100)")) + +par(mar=c(6, 5, 4, 2) + 0.1) + +sciplot::bargraph.CI(fattore, y, group = nomi.gruppi, data=dprof,legend=T,ylim=c(-2.4,1.6),x.leg=1.2, y.leg=-1.2,ylab="Punteggi standardizzati",xlab="",cex.lab=.8,cex.leg=.8, col=c("yellow","green","blue","red")) + +prop.table(table(ddef$gruppiattP, ddef$gruppiattM), margin = 1) +str(ddef$gruppiattP) + +cor(as.numeric(ddef$gruppiattP), as.numeric(ddef$gruppiattM)) + +#---- fit ext ---- + +hist(ddef$externalizing_sum) + +fit_ext <- glm(round(externalizing_sum,0) ~ genere * gruppiattP * gruppiattM, ddef, + family = "poisson") +summary(fit_ext) +car::Anova(fit_ext) +plot(fit_ext) + +performance::check_zeroinflation(fit_ext) + +plot(effects::allEffects(fit_ext)) + +# log transform +hist(log(ddef$externalizing_sum)) + +#---- fit ext no 0 ---- + +sum(ddef$externalizing_sum > 0) +ddef_ext0 <- ddef[ddef$externalizing_sum > 0, ] +hist(ddef_ext0$externalizing_sum) + +fit_ext0 <- glm(round(externalizing_sum,0) ~ genere * gruppiattP * gruppiattM, ddef_ext0, family = "poisson") +summary(fit_ext0) +car::Anova(fit_ext0) +plot(fit_ext0) + +plot(effects::allEffects(fit_ext0)) + + +summary(fit_ext0$residuals) +hist(fit_ext0$residuals) + + +#---- fit int ---- + +hist(ddef$internalizing_sum) +sum(ddef$internalizing_sum > 0) + +fit_int <- glm(round(internalizing_sum,0) ~ genere * gruppiattP * gruppiattM, ddef, family = "poisson") +summary(fit_int) +car::Anova(fit_int) + + +fit_int_ext <- lm(cbind(internalizing_sum, externalizing_sum) ~ genere * gruppiattP * gruppiattM, ddef) + + + +performance::check_zeroinflation(fit_int) + +plot(effects::allEffects(fit_int)) + + +#---- + + + diff --git a/Attachement.Rproj b/Attachement.Rproj new file mode 100644 index 0000000..88ff2b5 --- /dev/null +++ b/Attachement.Rproj @@ -0,0 +1,21 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: knitr +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageRoxygenize: rd,collate,namespace diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..7d58309 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,37 @@ +Package: Attachment +Type: Package +Title: What the Package Does (Title Case) +Version: 0.1.0 +Author@R: + c(person(given = "Claudio", + family = "Zandonella Callegher", + role = c("aut", "cre"), + email = "claudiozandonella@gmail.com", + comment = c(ORCID = "0000-0001-7721-6318")), + person(given = "Tatiana", + family = "Marci", + role = "aut", + email = "tatiana.marci@unipd.it", + comment = c(ORCID = "0000-0002-2813-0312")) + person(given = "Gianmarco", + family = "Altoè", + role = "aut", + email = "gianmarco.altoe@unipd.it", + comment = c(ORCID = "0000-0003-1154-9528")) + ) +Maintainer: Claudio Zandonella Callegher +Description: Analysis of atttachemnet theory +License: +Encoding: UTF-8 +Depends: + tidyverse, + devtools +Imports: + rmarkdown +Suggests: + knitr, + testthat (>= 3.0.0) +LazyData: true +Config/testthat/edition: 3 +RoxygenNote: 7.1.1 +VignetteBuilder: knitr diff --git a/Documents/Biblio-attachment.bib b/Documents/Biblio-attachment.bib new file mode 100644 index 0000000..52b59e5 --- /dev/null +++ b/Documents/Biblio-attachment.bib @@ -0,0 +1,408 @@ + +@article{bleiBayesianMixtureModels, + title = {Bayesian {{Mixture Models}} and the {{Gibbs Sampler}}}, + author = {Blei, David M}, + pages = {10}, + file = {/Users/claudio/Zotero/storage/HNGYXPA7/Blei - Bayesian Mixture Models and the Gibbs Sampler.pdf}, + langid = {english} +} + +@article{boing-messingBayesFactorsTesting2021, + title = {Bayes {{Factors}} for {{Testing Order Constraints}} on {{Variances}} of {{Dependent Outcomes}}}, + author = {Böing-Messing, Florian and Mulder, Joris}, + date = {2021-04-03}, + journaltitle = {The American Statistician}, + shortjournal = {The American Statistician}, + volume = {75}, + pages = {152--161}, + issn = {0003-1305, 1537-2731}, + doi = {10.1080/00031305.2020.1715257}, + url = {https://www.tandfonline.com/doi/full/10.1080/00031305.2020.1715257}, + urldate = {2021-05-11}, + abstract = {In statistical practice, researchers commonly focus on patterns in the means of multiple dependent outcomes while treating variances as nuisance parameters. However, in fact, there are often substantive reasons to expect certain patterns in the variances of dependent outcomes as well. For example, in a repeated measures study, one may expect the variance of the outcome to increase over time if the difference between subjects becomes more pronounced over time because the subjects respond differently to a given treatment. Such expectations can be formulated as order constrained hypotheses on the variances of the dependent outcomes. Currently, however, no methods exist for testing such hypotheses in a direct manner. To fill this gap, we develop a Bayes factor for this challenging testing problem. Our Bayes factor is based on the multivariate normal distribution with an unstructured covariance matrix, which is often used to model dependent outcomes. Order constrained hypotheses can then be formulated on the variances on the diagonal of the covariance matrix. To compute Bayes factors between multiple order constrained hypotheses, a prior distribution needs to be specified under every hypothesis to be tested. Here, we use the encompassing prior approach in which priors under order constrained hypotheses are truncations of the prior under the unconstrained hypothesis. The resulting Bayes factor is fully automatic in the sense that no subjective priors need to be specified by the user.}, + file = {/Users/claudio/Zotero/storage/4TI5Y359/Böing-Messing and Mulder - 2021 - Bayes Factors for Testing Order Constraints on Var.pdf}, + langid = {english}, + number = {2} +} + +@report{daganConfigurationsMotherChildFatherChild2021, + title = {Configurations of {{Mother}}-{{Child}} and {{Father}}-{{Child Attachment}} as {{Predictors}} of {{Internalizing}} and {{Externalizing Symptoms}}: {{An Individual Participant Data}} ({{IPD}}) {{Meta}}-{{Analysis}}}, + shorttitle = {Configurations of {{Mother}}-{{Child}} and {{Father}}-{{Child Attachment}} as {{Predictors}} of {{Internalizing}} and {{Externalizing Symptoms}}}, + author = {Dagan, Or and Schuengel, Carlo and Verhage, Marije and van IJzendoorn, Marinus H. and Sagi-Schwartz, Abraham and Madigan, Sheri and Duschinsky, Robbie and Roisman, Glenn I. and Bernard, Kristin and Bakermans-Kranenburg, Marian and Bureau, Jean-Francois and Volling, Brenda and Wong, Maria and Colonnesi, Cristina and Brown, Geoffrey and Eiden, Rina and Fearon, Pasco and Oosterman, Mirjam and Aviezer, Ora and Cummings, E. Mark}, + date = {2021-05-27}, + institution = {{PsyArXiv}}, + doi = {10.31234/osf.io/x4td2}, + url = {https://osf.io/x4td2}, + urldate = {2021-06-07}, + abstract = {An unsettled question in attachment theory and research is the extent to which children’s attachment patterns with mothers and fathers jointly predict developmental outcomes. In this study, we used individual participant data meta-analysis to assess whether early attachment networks with mothers and fathers are associated with children’s internalizing and externalizing symptoms. Following a pre-registered protocol, data from 9 studies and 1,097 children (mean age: 28.67 months) with attachment classifications to both mothers and fathers were included in analyses. We used a linear mixed effects analysis to assess differences in children’s internalizing and externalizing symptoms as assessed via the average of both maternal and paternal reports based on whether children had two, one, or no insecure (or disorganized) attachments. Results indicated that children with an insecure attachment relationship with one or both parents were at higher risk for elevated internalizing symptomatology compared with children who were securely attached to both parents. Children whose attachment relationships with both parents were classified as disorganized had more externalizing symptoms compared to children with either one or no disorganized attachment relationship with their parents. Across attachment classification networks and symptoms, findings suggest (a) a multiplicative effect when children have insecure or disorganized attachment to both parents, and (b) that mother-child and fatherchild attachment relationships may not differ in the roles they play in children’s development of internalizing and externalizing symptoms.}, + file = {/Users/claudio/Zotero/storage/MSZD47J9/Dagan et al. - 2021 - Configurations of Mother-Child and Father-Child At.pdf}, + langid = {english}, + options = {useprefix=true}, + type = {preprint} +} + +@online{DemystifyingJacobianAdjustment, + title = {Demystifying {{Jacobian}} Adjustment for Transformed Parameters in {{Stan}}}, + url = {http://rstudio-pubs-static.s3.amazonaws.com/486816_440106f76c944734a7d4c84761e37388.html}, + urldate = {2021-05-11}, + file = {/Users/claudio/Zotero/storage/KSCM9LRK/486816_440106f76c944734a7d4c84761e37388.html} +} + +@article{ghoshBayesianAnalysisZeroinflated2006, + title = {Bayesian Analysis of Zero-Inflated Regression Models}, + author = {Ghosh, Sujit K. and Mukhopadhyay, Pabak and Lu, Jye-Chyi(JC)}, + date = {2006-04}, + journaltitle = {Journal of Statistical Planning and Inference}, + shortjournal = {Journal of Statistical Planning and Inference}, + volume = {136}, + pages = {1360--1375}, + issn = {03783758}, + doi = {10.1016/j.jspi.2004.10.008}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0378375804004008}, + urldate = {2021-06-24}, + abstract = {In modeling defect counts collected from an established manufacturing processes, there are usually a relatively large number of zeros (non-defects). The commonly used models such as Poisson or Geometric distributions can underestimate the zero-defect probability and hence make it difficult to identify significant covariate effects to improve production quality. This article introduces a flexible class of zero inflated models which includes other familiar models such as the Zero Inflated Poisson (ZIP) models, as special cases. A Bayesian estimation method is developed as an alternative to traditionally used maximum likelihood based methods to analyze such data. Simulation studies show that the proposed method has better finite sample performance than the classical method with tighter interval estimates and better coverage probabilities. A real-life data set is analyzed to illustrate the practicability of the proposed method easily implemented using WinBUGS.}, + file = {/Users/claudio/Zotero/storage/GKHC66BJ/Ghosh et al. - 2006 - Bayesian analysis of zero-inflated regression mode.pdf}, + langid = {english}, + number = {4} +} + +@online{gronauBridgesamplingPackageEstimating2018, + title = {Bridgesampling: {{An R Package}} for {{Estimating Normalizing Constants}}}, + shorttitle = {Bridgesampling}, + author = {Gronau, Quentin F. and Singmann, Henrik and Wagenmakers, Eric-Jan}, + date = {2018-09-18}, + url = {http://arxiv.org/abs/1710.08162}, + urldate = {2021-05-27}, + abstract = {Statistical procedures such as Bayes factor model selection and Bayesian model averaging require the computation of normalizing constants (e.g., marginal likelihoods). These normalizing constants are notoriously difficult to obtain, as they usually involve highdimensional integrals that cannot be solved analytically. Here we introduce an R package that uses bridge sampling (Meng and Wong 1996; Meng and Schilling 2002) to estimate normalizing constants in a generic and easy-to-use fashion. For models implemented in Stan, the estimation procedure is automatic. We illustrate the functionality of the package with three examples.}, + archiveprefix = {arXiv}, + eprint = {1710.08162}, + eprinttype = {arxiv}, + file = {/Users/claudio/Zotero/storage/8JED3CYZ/Gronau et al. - 2018 - bridgesampling An R Package for Estimating Normal.pdf}, + keywords = {Statistics - Computation}, + langid = {english}, + primaryclass = {stat} +} + +@article{gronauTutorialBridgeSampling2017, + title = {A Tutorial on Bridge Sampling}, + author = {Gronau, Quentin F. and Sarafoglou, Alexandra and Matzke, Dora and Ly, Alexander and Boehm, Udo and Marsman, Maarten and Leslie, David S. and Forster, Jonathan J. and Wagenmakers, Eric-Jan and Steingroever, Helen}, + date = {2017-12}, + journaltitle = {Journal of Mathematical Psychology}, + shortjournal = {Journal of Mathematical Psychology}, + volume = {81}, + pages = {80--97}, + issn = {00222496}, + doi = {10.1016/j.jmp.2017.09.005}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0022249617300640}, + urldate = {2021-05-27}, + abstract = {The marginal likelihood plays an important role in many areas of Bayesian statistics such as parameter estimation, model comparison, and model averaging. In most applications, however, the marginal likelihood is not analytically tractable and must be approximated using numerical methods. Here we provide a tutorial on bridge sampling (Bennett, 1976; Meng \& Wong, 1996), a reliable and relatively straightforward sampling method that allows researchers to obtain the marginal likelihood for models of varying complexity. First, we introduce bridge sampling and three related sampling methods using the beta-binomial model as a running example. We then apply bridge sampling to estimate the marginal likelihood for the Expectancy Valence (EV) model—a popular model for reinforcement learning. Our results indicate that bridge sampling provides accurate estimates for both a single participant and a hierarchical version of the EV model. We conclude that bridge sampling is an attractive method for mathematical psychologists who typically aim to approximate the marginal likelihood for a limited set of possibly high-dimensional models.}, + file = {/Users/claudio/Zotero/storage/ZREL5TQV/Gronau et al. - 2017 - A tutorial on bridge sampling.pdf}, + langid = {english} +} + +@article{guApproximatedAdjustedFractional2018, + title = {Approximated Adjusted Fractional {{Bayes}} Factors: {{A}} General Method for Testing Informative Hypotheses}, + shorttitle = {Approximated Adjusted Fractional {{Bayes}} Factors}, + author = {Gu, Xin and Mulder, Joris and Hoijtink, Herbert}, + date = {2018-05}, + journaltitle = {British Journal of Mathematical and Statistical Psychology}, + shortjournal = {Br J Math Stat Psychol}, + volume = {71}, + pages = {229--261}, + issn = {00071102}, + doi = {10.1111/bmsp.12110}, + url = {http://doi.wiley.com/10.1111/bmsp.12110}, + urldate = {2021-05-13}, + file = {/Users/claudio/Zotero/storage/X623DBMS/Gu et al. - 2018 - Approximated adjusted fractional Bayes factors A .pdf}, + keywords = {Letto}, + langid = {english}, + number = {2} +} + +@article{guBayesianEvaluationInequality2014, + title = {Bayesian Evaluation of Inequality Constrained Hypotheses.}, + author = {Gu, Xin and Mulder, Joris and Deković, Maja and Hoijtink, Herbert}, + date = {2014}, + journaltitle = {Psychological Methods}, + shortjournal = {Psychological Methods}, + volume = {19}, + pages = {511--527}, + issn = {1939-1463, 1082-989X}, + doi = {10.1037/met0000017}, + url = {http://doi.apa.org/getdoi.cfm?doi=10.1037/met0000017}, + urldate = {2021-05-13}, + abstract = {Bayesian evaluation of inequality constrained hypotheses enables researchers to investigate their expectations with respect to the structure among model parameters. This article proposes an approximate Bayes procedure that can be used for the selection of the best of a set of inequality constrained hypotheses based on the Bayes factor in a very general class of statistical models. The software package BIG is provided such that psychologists can use the approach proposed for the analysis of their own data. To illustrate the approximate Bayes procedure and the use of BIG, we evaluate inequality constrained hypotheses in a path model and a logistic regression model. Two simulation studies on the performance of our approximate Bayes procedure show that it results in accurate Bayes factors.}, + file = {/Users/claudio/Zotero/storage/PUHJAVZD/Gu et al. - 2014 - Bayesian evaluation of inequality constrained hypo.pdf}, + keywords = {Letto}, + langid = {english}, + number = {4} +} + +@article{heckMultinomialModelsLinear2019, + title = {Multinomial Models with Linear Inequality Constraints: {{Overview}} and Improvements of Computational Methods for {{Bayesian}} Inference}, + shorttitle = {Multinomial Models with Linear Inequality Constraints}, + author = {Heck, Daniel W. and Davis-Stober, Clintin P.}, + date = {2019-08}, + journaltitle = {Journal of Mathematical Psychology}, + shortjournal = {Journal of Mathematical Psychology}, + volume = {91}, + pages = {70--87}, + issn = {00222496}, + doi = {10.1016/j.jmp.2019.03.004}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0022249618301457}, + urldate = {2021-05-11}, + abstract = {Many psychological theories can be operationalized as linear inequality constraints on the parameters of multinomial distributions (e.g., discrete choice analysis). These constraints can be described in two equivalent ways: Either as the solution set to a system of linear inequalities or as the convex hull of a set of extremal points (vertices). For both representations, we describe a general Gibbs sampler for drawing posterior samples in order to carry out Bayesian analyses. We also summarize alternative sampling methods for estimating Bayes factors for these model representations using the encompassing Bayes factor method. We introduce the R package multinomineq, which provides an easily-accessible interface to a computationally efficient implementation of these techniques.}, + file = {/Users/claudio/MEGA/Zotero/Heck_Davis-Stober_2019_Multinomial models with linear inequality constraints.pdf}, + keywords = {To read}, + langid = {english} +} + +@report{heckReviewApplicationsBayes2020, + title = {A {{Review}} of {{Applications}} of the {{Bayes Factor}} in {{Psychological Research}}}, + author = {Heck, Daniel W. and Boehm, Udo and Böing-Messing, Florian and Bürkner, Paul - Christian and Derks, Koen and Dienes, Zoltan and Fu, Qianrao and Gu, Xin and Karimova, Diana and Kiers, Henk and Klugkist, Irene and Kuiper, Rebecca M. and Lee, Michael David and Leenders, Roger and Leplaa, Hidde Jelmer and Linde, Maximilian and Ly, Alexander and Meijerink-Bosman, Marlyne and Moerbeek, Mirjam and Mulder, Joris and Palfi, Bence and Schönbrodt, Felix D. and Tendeiro, Jorge and van den Bergh, Don and Van Lissa, Caspar J. and van Ravenzwaaij, Don and {vanpaemel}, wolf and Wagenmakers, Eric-Jan and Williams, Donald Ray and Zondervan-Zwijnenburg, Marielle and Hoijtink, Herbert}, + date = {2020-10-20}, + institution = {{PsyArXiv}}, + doi = {10.31234/osf.io/cu43g}, + url = {https://osf.io/cu43g}, + urldate = {2021-05-31}, + abstract = {The last 25 years have shown a steady increase in attention for the Bayes factor as a tool for hypothesis evaluation and model selection. The present review highlights the potential of the Bayes factor in psychological research. We discuss six types of applications: Bayesian evaluation of point null, interval, and informative hypotheses, Bayesian evidence synthesis, Bayesian variable selection and model averaging, and Bayesian evaluation of cognitive models. We elaborate what each application entails, give illustrative examples, and provide an overview of key references and software with links to other applications. The paper is concluded with a discussion of the opportunities and pitfalls of Bayes factor applications and a sketch of corresponding future research lines.}, + file = {/Users/claudio/Zotero/storage/FTL4BS3W/Heck et al. - 2020 - A Review of Applications of the Bayes Factor in Ps.pdf}, + langid = {english}, + options = {useprefix=true}, + type = {preprint} +} + +@book{hoijtinkInformativeHypotheses2011, + title = {Informative {{Hypotheses}}}, + author = {Hoijtink, Herbert}, + date = {2011-10-26}, + edition = {0}, + publisher = {{Chapman and Hall/CRC}}, + doi = {10.1201/b11158}, + url = {https://www.taylorfrancis.com/books/9781439880524}, + urldate = {2021-05-31}, + file = {/Users/claudio/Zotero/storage/VQ4YJNCD/Hoijtink - 2011 - Informative Hypotheses.pdf}, + isbn = {978-1-4398-8052-4}, + langid = {english} +} + +@article{hoijtinkTutorialTestingHypotheses2019, + title = {A Tutorial on Testing Hypotheses Using the {{Bayes}} Factor.}, + author = {Hoijtink, Herbert and Mulder, Joris and van Lissa, Caspar and Gu, Xin}, + date = {2019-10}, + journaltitle = {Psychological Methods}, + shortjournal = {Psychological Methods}, + volume = {24}, + pages = {539--556}, + issn = {1939-1463, 1082-989X}, + doi = {10.1037/met0000201}, + url = {http://doi.apa.org/getdoi.cfm?doi=10.1037/met0000201}, + urldate = {2021-05-13}, + abstract = {Learning about hypothesis evaluation using the Bayes factor could enhance psychological research. In contrast to null-hypothesis significance testing: it renders the evidence in favor of each of the hypotheses under consideration (it can be used to quantify support for the null-hypothesis) instead of a dichotomous reject/do-not-reject decision; it can straightforwardly be used for the evaluation of multiple hypotheses without having to bother about the proper manner to account for multiple testing; and, it allows continuous re-evaluation of hypotheses after additional data have been collected (Bayesian updating). This tutorial addresses researchers considering to evaluate their hypotheses by means of the Bayes factor. The focus is completely applied and each topic discussed is illustrated using Bayes factors for the evaluation of hypotheses in the context of an ANOVA model, obtained using the R package bain. Readers can execute all the analyses presented while reading this tutorial if they download bain and the R-codes used. It will be elaborated in a completely non-technical manner: what the Bayes factor is, how it can be obtained, how Bayes factors should be interpreted, and what can be done with Bayes factors. After reading this tutorial and executing the associated code, researchers will be able to use their own data for the evaluation of hypotheses by means of the Bayes factor, not only in the context of ANOVA models, but also in the context of other statistical models.}, + file = {/Users/claudio/Zotero/storage/UH74J8BX/Hoijtink et al. - 2019 - A tutorial on testing hypotheses using the Bayes f.pdf}, + langid = {english}, + number = {5}, + options = {useprefix=true} +} + +@article{katoBayesianApproachInequality2006, + title = {A {{Bayesian}} Approach to Inequality Constrained Linear Mixed Models: Estimation and Model Selection}, + shorttitle = {A {{Bayesian}} Approach to Inequality Constrained Linear Mixed Models}, + author = {Kato, Bernet S and Hoijtink, Herbert}, + date = {2006-10}, + journaltitle = {Statistical Modelling}, + shortjournal = {Statistical Modelling}, + volume = {6}, + pages = {231--249}, + issn = {1471-082X, 1477-0342}, + doi = {10.1191/1471082X06st119oa}, + url = {http://journals.sagepub.com/doi/10.1191/1471082X06st119oa}, + urldate = {2021-05-11}, + abstract = {Constrained parameter problems arise in a wide variety of applications. This article deals with estimation and model selection in linear mixed models with inequality constraints on the parameters. It is shown that different theories can be translated into statistical models by putting constraints on the model parameters yielding a set of competing models. A new approach based on the principle of encompassing priors is proposed and used to compute Bayes factors and subsequently posterior model probabilities. Model selection is based on posterior model probabilities. The approach is illustrated using a longitudinal data set.}, + file = {/Users/claudio/MEGA/Zotero/Kato_Hoijtink_2006_A Bayesian approach to inequality constrained linear mixed models.pdf}, + keywords = {Letto}, + langid = {english}, + number = {3} +} + +@article{liuBAYESIANINFERENCEZEROINFLATED, + title = {{{BAYESIAN INFERENCE FOR ZERO}}-{{INFLATED POISSON REGRESSION MODELS}}}, + author = {Liu, Hui and Powers, Daniel A}, + pages = {34}, + abstract = {Count data with excess zeros are common in social science research and can be considered as a special case of mixture structured data. We exploit the flexibility of the Bayesian analytic approach to model the mixture data structure inherent in zero-inflated count data by using the zero-inflated Poisson (ZIP) model. We discuss the importance of modelling excess-zero count data in social sciences and review the distributional properties of zero-inflated count data, with special attention given to its mixture data structure in the context of Bayesian modelling. We illustrate the methodology using data from the Americans’ changing lives (ACL) survey on cigarette smoking. Results from predictive checks suggest that the proposed Bayesian ZIP model provides a good fit to the 2010 Mathematics Subject Classification: 62F15, 62J12.}, + file = {/Users/claudio/Zotero/storage/IJ9BAQ6U/Liu and Powers - BAYESIAN INFERENCE FOR ZERO-INFLATED POISSON REGRE.pdf}, + keywords = {To read}, + langid = {english} +} + +@article{liuZeroandoneinflatedPoissonRegression2021, + title = {Zero-and-One-Inflated {{Poisson}} Regression Model}, + author = {Liu, Wenchen and Tang, Yincai and Xu, Ancha}, + date = {2021-04}, + journaltitle = {Statistical Papers}, + shortjournal = {Stat Papers}, + volume = {62}, + pages = {915--934}, + issn = {0932-5026, 1613-9798}, + doi = {10.1007/s00362-019-01118-7}, + url = {http://link.springer.com/10.1007/s00362-019-01118-7}, + urldate = {2021-06-24}, + abstract = {In this paper, a zero-and-one-inflated Poisson (ZOIP) regression model is proposed. The maximum likelihood estimation (MLE) and Bayesian estimation for this model are investigated. Three estimation methods of the ZOIP regression model are obtained based on data augmentation method which is expectation-maximization (EM) algorithm, generalized expectation-maximization (GEM) algorithm and Gibbs sampling respectively. A simulation study is conducted to assess the performance of the proposed estimation for various sample sizes. Finally, an accidental deaths data set is analyzed to illustrate the practicability of the proposed method.}, + file = {/Users/claudio/Zotero/storage/QFQFYRN8/Liu et al. - 2021 - Zero-and-one-inflated Poisson regression model.pdf}, + langid = {english}, + number = {2} +} + +@article{moreyPhilosophyBayesFactors2016, + title = {The Philosophy of {{Bayes}} Factors and the Quantification of Statistical Evidence}, + author = {Morey, Richard D. and Romeijn, Jan-Willem and Rouder, Jeffrey N.}, + date = {2016-06}, + journaltitle = {Journal of Mathematical Psychology}, + shortjournal = {Journal of Mathematical Psychology}, + volume = {72}, + pages = {6--18}, + issn = {00222496}, + doi = {10.1016/j.jmp.2015.11.001}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0022249615000723}, + urldate = {2021-05-13}, + abstract = {A core aspect of science is using data to assess the degree to which data provide evidence for competing claims, hypotheses, or theories. Evidence is by definition something that should change the credibility of a claim in a reasonable person’s mind. However, common statistics, such as significance testing and confidence intervals have no interface with concepts of belief, and thus it is unclear how they relate to statistical evidence. We explore the concept of statistical evidence, and how it can be quantified using the Bayes factor. We also discuss the philosophical issues inherent in the use of the Bayes factor.}, + file = {/Users/claudio/Zotero/storage/35H3EWNI/Morey et al. - 2016 - The philosophy of Bayes factors and the quantifica.pdf}, + langid = {english} +} + +@article{mulderBayesFactorsTesting2016, + title = {Bayes Factors for Testing Order-Constrained Hypotheses on Correlations}, + author = {Mulder, Joris}, + date = {2016-06}, + journaltitle = {Journal of Mathematical Psychology}, + shortjournal = {Journal of Mathematical Psychology}, + volume = {72}, + pages = {104--115}, + issn = {00222496}, + doi = {10.1016/j.jmp.2014.09.004}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0022249614000601}, + urldate = {2021-05-13}, + abstract = {Correlation coefficients play a key role in the social and behavioral sciences for quantifying the degree of linear association between variables. A Bayes factor is proposed that allows researchers to test hypotheses with order constraints on correlation coefficients in a direct manner. This Bayes factor balances between fit and complexity of order-constrained hypotheses in a natural way. A diffuse prior on the correlation matrix is used that minimizes prior shrinkage and results in most evidence for an order-constrained hypothesis that is supported by the data. An efficient method is proposed for the computation of the Bayes factor. A key aspect in the computation is a Fisher Z transformation on the posterior distribution of the correlations such that an approximately normal distribution is obtained. The methodology is implemented in a freely downloadable software program called ‘‘BOCOR’’. The methods are applied to a multitrait–multimethod analysis, a repeated measures study, and a study on directed moderator effects. © 2014 Elsevier Inc. All rights reserved.}, + file = {/Users/claudio/Zotero/storage/PYUHF2IL/Mulder - 2016 - Bayes factors for testing order-constrained hypoth.pdf}, + keywords = {Letto Parzialmente}, + langid = {english} +} + +@online{mulderBayesFactorTesting2019, + title = {Bayes Factor Testing of Equality and Order Constraints on Measures of Association in Social Research}, + author = {Mulder, Joris and Gelissen, John P. T. M.}, + date = {2019-04-03}, + url = {http://arxiv.org/abs/1807.05819}, + urldate = {2021-05-11}, + abstract = {Measures of association play a central role in the social sciences to quantify the strength of a linear relationship between the variables of interest. In many applications researchers can translate scientific expectations to hypotheses with equality and/or order constraints on these measures of association. In this paper a Bayes factor test is proposed for testing multiple hypotheses with constraints on the measures of association between ordinal and/or continuous variables, possibly after correcting for certain covariates. This test can be used to obtain a direct answer to the research question how much evidence there is in the data for a social science theory relative to competing theories. The accompanying software package ‘BCT’ allows users to apply the methodology in an easy manner. An empirical application from leisure studies about the associations between life, leisure and relationship satisfaction and an application about the differences about egalitarian justice beliefs across countries are used to illustrate the methodology.}, + archiveprefix = {arXiv}, + eprint = {1807.05819}, + eprinttype = {arxiv}, + file = {/Users/claudio/Zotero/storage/D2X73UGU/Mulder and Gelissen - 2019 - Bayes factor testing of equality and order constra.pdf}, + keywords = {Letto Parzialmente,Statistics - Methodology}, + langid = {english}, + primaryclass = {stat} +} + +@article{mulderSimpleBayesianTesting2019, + title = {Simple {{Bayesian}} Testing of Scientific Expectations in Linear Regression Models}, + author = {Mulder, J. and Olsson-Collentine, A.}, + date = {2019-06}, + journaltitle = {Behavior Research Methods}, + shortjournal = {Behav Res}, + volume = {51}, + pages = {1117--1130}, + issn = {1554-3528}, + doi = {10.3758/s13428-018-01196-9}, + url = {http://link.springer.com/10.3758/s13428-018-01196-9}, + urldate = {2021-06-07}, + abstract = {Scientific theories can often be formulated using equality and order constraints on the relative effects in a linear regression model. For example, it may be expected that the effect of the first predictor is larger than the effect of the second predictor, and the second predictor is expected to be larger than the third predictor. The goal is then to test such expectations against competing scientific expectations or theories. In this paper, a simple default Bayes factor test is proposed for testing multiple hypotheses with equality and order constraints on the effects of interest. The proposed testing criterion can be computed without requiring external prior information about the expected effects before observing the data. The method is implemented in R-package called ‘lmhyp’ which is freely downloadable and ready to use. The usability of the method and software is illustrated using empirical applications from the social and behavioral sciences.}, + file = {/Users/claudio/MEGA/Zotero/Mulder_Olsson-Collentine_2019_Simple Bayesian testing of scientific expectations in linear regression models.pdf}, + langid = {english}, + number = {3} +} + +@article{rodrigues-mottaZeroInflatedPoissonModel2007, + title = {A {{Zero}}-{{Inflated Poisson Model}} for {{Genetic Analysis}} of the {{Number}} of {{Mastitis Cases}} in {{Norwegian Red Cows}}}, + author = {Rodrigues-Motta, M. and Gianola, D. and Heringstad, B. and Rosa, G.J.M. and Chang, Y.M.}, + date = {2007-11}, + journaltitle = {Journal of Dairy Science}, + shortjournal = {Journal of Dairy Science}, + volume = {90}, + pages = {5306--5315}, + issn = {00220302}, + doi = {10.3168/jds.2006-898}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0022030207720003}, + urldate = {2021-05-31}, + abstract = {The objective was to extend a zero-inflated Poisson (ZIP) model to account for correlated genetic effects, and to use this model to analyze the number of clinical mastitis cases in Norwegian Red cows. The ZIP model is suitable for analysis of count data containing an excess of zeros relative to what is expected from Poisson sampling. A ZIP model was developed and compared with a corresponding Poisson model. The Poisson parameter followed a hierarchical structure, and a residual term accounting for overdispersion was included. In both models, the Poisson parameter was regressed 1) on the year, month, and age at first calving; 2) on the logarithm of the number of days elapsed from calving to the end of first lactation; and c) on herd and sire effects. Herd and sire effects were assigned normal prior distributions in a Bayesian analysis, corresponding to a random effects treatment in a frequentist analysis. An analysis of residuals favored the Poisson model when there were 2 or more cases of mastitis during first lactation, with very small differences between the ZIP and Poisson models at 0 and 1 cases. However, the residual assessment was not satisfactory for either of the models. The ZIP model, on the other hand, had a better predictive ability than the corresponding Poisson model. Posterior means of the sire, herd, and residual variances in the ZIP model (log scale) were 0.09, 0.37, and 0.36, respectively, highlighting the importance of herds as a source of variation in clinical mastitis. The correlation between sire rankings from the ZIP and Poisson models was 0.98. A weaker correlation would be expected in a population with more severe inflation at zero than the present one. The estimate of the perfect state probability p was 0.32, indicating that 32\% of the animals would be in the perfect state, either because they are resistant or because they were not exposed to mastitis.}, + file = {/Users/claudio/Zotero/storage/3YUADYFB/Rodrigues-Motta et al. - 2007 - A Zero-Inflated Poisson Model for Genetic Analysis.pdf}, + keywords = {Letto Parzialmente}, + langid = {english}, + number = {11} +} + +@article{rodriguesBayesianAnalysisZeroInflated2003, + title = {Bayesian {{Analysis}} of {{Zero}}-{{Inflated Distributions}}}, + author = {Rodrigues, Josemar}, + date = {2003-01-03}, + journaltitle = {Communications in Statistics - Theory and Methods}, + shortjournal = {Communications in Statistics - Theory and Methods}, + volume = {32}, + pages = {281--289}, + issn = {0361-0926, 1532-415X}, + doi = {10.1081/STA-120018186}, + url = {http://www.tandfonline.com/doi/abs/10.1081/STA-120018186}, + urldate = {2021-05-28}, + abstract = {In this paper zero-inflated distributions (ZID) are studied from the Bayesian point of view using the data augmentation algorithm. This type of discrete model arises in count data with excess of zeros. The zero-inflated Poisson distribution (ZIP) and an illustrative example via MCMC algorithm are considered.}, + file = {/Users/claudio/Zotero/storage/YFDYQKLH/Rodrigues - 2003 - Bayesian Analysis of Zero-Inflated Distributions.pdf}, + langid = {english}, + number = {2} +} + +@article{rodriguez-yamEfficientGibbsSampling, + title = {Efficient {{Gibbs Sampling}} for {{Constrained Linear Regression}}}, + author = {Rodriguez-Yam, Gabriel and Davis, Richard A and Scharf, Louis L}, + pages = {6}, + abstract = {In this paper we consider parameter estimation in a linear regression setting with inequality linear constraints on the regression parameters. Most other research on this topic has typically been addressed from a Bayesian perspective. Geweke (1996) and Chen and Deely (1996) implemented the Gibbs sampler to generate samples from the posterior distribution. However, these implementations can often exhibit poor mixing and slow convergence. This paper overcomes these limitations with a new implementation of the Gibbs sampler. In addition, this procedure allows for the number of constraints to exceed the parameter dimension and is able to cope with equality linear constraints.}, + file = {/Users/claudio/Zotero/storage/XTSV2X3F/Rodriguez-Yam et al. - Efficient Gibbs Sampling for Constrained Linear Regr.pdf}, + langid = {english} +} + +@article{vandeschootIntroductionBayesianModel2011, + title = {An Introduction to {{Bayesian}} Model Selection for Evaluating Informative Hypotheses}, + author = {van de Schoot, Rens and Mulder, Joris and Hoijtink, Herbert and Van Aken, Marcel A. G. and Semon Dubas, Judith and Orobio de Castro, Bram and Meeus, Wim and Romeijn, Jan-Willem}, + date = {2011-11}, + journaltitle = {European Journal of Developmental Psychology}, + volume = {8}, + pages = {713--729}, + issn = {1740-5629, 1740-5610}, + doi = {10.1080/17405629.2011.621799}, + url = {http://www.tandfonline.com/doi/abs/10.1080/17405629.2011.621799}, + urldate = {2019-05-22}, + file = {/Users/claudio/MEGA/Zotero/van de Schoot et al_2011_An introduction to Bayesian model selection for evaluating informative.pdf;/Users/claudio/MEGA/Zotero/van de Schoot et al_2011_An introduction to Bayesian model selection for evaluating informative2.pdf}, + keywords = {Bayesian,Informative Hypotesis,Letto,Model Comparison}, + langid = {english}, + number = {6}, + options = {useprefix=true} +} + +@article{wagenmakersBayesianHypothesisTesting2010, + title = {Bayesian Hypothesis Testing for Psychologists: {{A}} Tutorial on the {{Savage}}–{{Dickey}} Method}, + shorttitle = {Bayesian Hypothesis Testing for Psychologists}, + author = {Wagenmakers, Eric-Jan and Lodewyckx, Tom and Kuriyal, Himanshu and Grasman, Raoul}, + date = {2010-05}, + journaltitle = {Cognitive Psychology}, + shortjournal = {Cognitive Psychology}, + volume = {60}, + pages = {158--189}, + issn = {00100285}, + doi = {10.1016/j.cogpsych.2009.12.001}, + url = {https://linkinghub.elsevier.com/retrieve/pii/S0010028509000826}, + urldate = {2021-05-25}, + abstract = {In the field of cognitive psychology, the p-value hypothesis test has established a stranglehold on statistical reporting. This is unfortunate, as the p-value provides at best a rough estimate of the evidence that the data provide for the presence of an experimental effect. An alternative and arguably more appropriate measure of evidence is conveyed by a Bayesian hypothesis test, which prefers the model with the highest average likelihood. One of the main problems with this Bayesian hypothesis test, however, is that it often requires relatively sophisticated numerical methods for its computation. Here we draw attention to the Savage–Dickey density ratio method, a method that can be used to compute the result of a Bayesian hypothesis test for nested models and under certain plausible restrictions on the parameter priors. Practical examples demonstrate the method’s validity, generality, and flexibility.}, + file = {/Users/claudio/Zotero/storage/YDPVNINF/Wagenmakers et al. - 2010 - Bayesian hypothesis testing for psychologists A t.pdf}, + langid = {english}, + number = {3} +} + + diff --git a/Documents/Report-Project/Attachment_Analysis_Report.Rmd b/Documents/Report-Project/Attachment_Analysis_Report.Rmd new file mode 100644 index 0000000..eab378b --- /dev/null +++ b/Documents/Report-Project/Attachment_Analysis_Report.Rmd @@ -0,0 +1,458 @@ +--- +title: 'Attachment Analysis: Report Project' +author: 'CZC, TM & GA' +date: "`r format(Sys.time(), '%d %B, %Y')`" +output: + bookdown::html_document2: + toc: true + toc_depth: 3 + toc_float: true + collapsed: false + css: ["css/custom.css"] +linestretch: 1.5 +number_sections: true +bibliography: "../Biblio-attachment.bib" +csl: apa.csl +link-citations: true +editor_options: + chunk_output_type: console +--- + + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = FALSE, collapse = TRUE) + +library(tidyverse) +library(kableExtra) +library(mclust) + +devtools::load_all("../../") + +path_drake <- "../../.drake" + +#---- load drake ---- + +# Data +drake::loadd(data_munged) + +# Cluster +drake::loadd(cluster_mother_tat) +drake::loadd(cluster_father_tat) +drake::loadd(data_cluster_tat) +drake::loadd(mclust_mother) +drake::loadd(mclust_father) + +# brms models +drake::loadd(fit_int) +``` + + +# Descriptive Statistics + + +## Dependent Variable + +### Externalizing +```{r} +ggplot(data_munged) + + geom_bar(aes(x = externalizing_sum, fill = genere), alpha = 0.7) + + facet_grid(genere ~ .) +``` + +### Internalizing +```{r} +ggplot(data_munged) + + geom_bar(aes(x = internalizing_sum, fill = genere), alpha = 0.7) + + facet_grid(genere ~ .) +``` + + +# Cluster Analysis + +## Mother + +```{r} +plot(cluster_mother_tat, main = "Dendrogramma") +rect.hclust(cluster_mother_tat, k=4, border="red") +``` + +Frequencies + +```{r, echo = TRUE} +table(data_cluster_tat$cluster_mother_tat) + +prop.table(table(data_cluster_tat$cluster_mother_tat)) + +prop.table(table(data_cluster_tat$cluster_mother_tat, + data_cluster_tat$genere),2) + +summary(table(data_cluster_tat$cluster_mother_tat, + data_cluster_tat$genere)) +``` + +Plot Scores + +```{r} +plot_scores_mother(data = data_cluster_tat) +``` + +#### Mclust check + +```{r} +plot(mclust_mother) +summary(mclust_mother) +``` + + +## Father + +```{r} +plot(cluster_father_tat, main = "Dendrogramma") +rect.hclust(cluster_father_tat, k=4, border="red") +``` + +Frequencies + +```{r, echo = TRUE} +table(data_cluster_tat$cluster_father_tat) + +prop.table(table(data_cluster_tat$cluster_father_tat)) + +prop.table(table(data_cluster_tat$cluster_father_tat, + data_cluster_tat$genere),2) + +summary(table(data_cluster_tat$cluster_father_tat, + data_cluster_tat$genere)) +``` + +Plot Scores + +```{r} +plot_scores_father(data = data_cluster_tat) +``` + +#### Mclust check + +```{r} +plot(mclust_father) +summary(mclust_father) +``` + +# Statistical Approach + +In this section we describe the statistical approach. First we define the statistical model used to describe the data generative process. Subsequently, the approach for testing hypotheses with inequality and equality constraints is discussed. Finally, the Gibbs sampler to estimate the model with inequality constraints is presented. + +## Model Definition + +### Psychologial Background + +In this analysis, we want to evaluate children psychological problems according to children gender and children-mother and children-father attachment relationship. Note that attachment relationship is considered as a categorical variable. + +### Zero Inflated Poisson (ZIP) {#model} + + +The dependent variable is a positive discrete variable. Considering Figure \@ref(fig:plot-int), it is possible to observe the the distribution of the dependent variable is characterized by a zero inflation. + +```{r plot-int, fig.cap = "Dependent variable distribution"} +ggplot(data_munged) + + geom_bar(aes(x = internalizing_sum, fill = 1), alpha = 0.7) + + labs(x = "Internalizing Problems") + + theme_bw() + + theme(legend.position = "none") +``` + +Thus, we decided to model the data considering a Zero Inflated Poisson (ZIP): + +$$ +y_i \sim ZIPoisson(p_i, \lambda_i) +$$ + +Where, $1 - p_i$ indicates the probability of given observation $y_i$ being generated from a Poisson distribution with mean $\lambda_i$. In particular, + +$$ +p_i = \text{logit}^{-1}(X_i^T\beta_{p}),\\ +\lambda_i = \text{exp}(X_i^T\beta_{\lambda}). +$$ + +That is, we want to model both $p_i$ and $\lambda_i$ according to the variables previously defined estimating respectively parameters vector $\beta_p$ and $\beta_{\lambda}$. Note that these are the only parameters in the model. + +This model can be easily estimate using `brms`. However, currently it is not possible to include inequality and equality constraints in the parameters prior definition. + +Note that, the actual distribution of the dependent variable is truncated (limit given by the maximum questionnaire score). However, we do not consider this aspect to avoid a further level of complexity in the model. + +## Test of Hypotheses with Inequality and Equality Constraints + +**Bayes Factor** to evaluate different hypotheses with inequality and equality constraints can be obtained following the **Encompassing Prior Approach**. That is, an unconstrained model is considered as reference with encompassing priors that are uninformative with respect to the parameters. All other models are nested in this model and the priors ore obtained by restricting the parameter space in accordance with the constraints imposed by each model. + +### Inequality Constraints: A First Approach + +@katoBayesianApproachInequality2006 presented an approach to evaluate **inequality constraints** based on the principle of encompassing priors. Given an unconstrained model $M_u$ and a nested model $M_i$, the prior and posterior densities of $M_i$ can be rewritten in terms of the prior and posterior densities of $M_u$: + +$$ +\pi(\theta|M_i) = B_i \times \pi(\theta|M_u),\\ +Pr(\theta|D,\ M_i) = A_i \times Pr(\theta|D,\ M_u). +$$ + +Where, $1/B_i$ indicates the proportion of the prior distribution of $M_u$ in agreement with $M_i$ and $1/A_i$ is the proportion of the posterior distribution of $M_u$ in agreement with $M_i$. + +Assuming equal prior probability of the two models ($\pi(M_u) = \pi(M_i)$), the Bayes Factor is obtained as follow: + +\begin{aligned} +BF_{iu} =& \frac{Pr(D|M_i)}{Pr(D|M_u)} =\\ \\ +&\frac{L(D|\theta)\pi(\theta|M_i)/Pr(\theta|D,\ M_i)}{L(D|\theta)\pi(\theta|M_u)/Pr(\theta|D,\ M_u)} = \\ \\ +&\frac{\pi(\theta|M_i)Pr(\theta|D,\ M_u)}{\pi(\theta|M_u)Pr(\theta|D,\ M_i)} = \\ \\ +&\frac{[B_i \times \pi(\theta|M_u)] Pr(\theta|D,\ M_u)}{\pi(\theta|M_u)[A_i \times Pr(\theta|D,\ M_u)]} = \frac{B_i}{A_i} +\end{aligned} + +That is, the Bayes Factor is the ratio between $1/A_i$ (i.e. the proportion of the posterior distribution of $M_u$ in agreement with $M_i$ that provides information about the *fit* of $M_i$ relative to $M_u$), and $1/B_i$ (i.e., the proportion of the prior distribution of $M_u$ in agreement with $M_i$ that provides information about the *complexity* of $M_i$ relative to $M_u$). + +Note that Bayes Factor is obtained without computing the marginal likelihoods of the two models, but only considering the proportions of the prior and the posterior in agreement with constraints of each of the models. + +#### Pro {-} + +- Given the prior and posterior samples of the model, it is straightforward to compute $1/A_i$ and $1/B_i$. The analytical distributions are not required. +- If Prior distributions of $\beta_p$ and $\beta_{\lambda}$ are vague enough and sample is large, influence of the choice of the encompassing prior on the posterior distributions is negligible, and subsequently also its influence on $1/A_i$ is negligible. +- For constraints of the form $\beta_i > \beta_j > \beta_k$ with $\beta_{i,j,k} \stackrel{i.i.d.}{\sim} \pi(\beta)$, the prior probability of the constraints is independent of the actual choice of the encompassing prior, and subsequently also $1/B_i$ is not influenced. + +#### Cons {-} + +- It is not possible to evaluate **equality constraints**. In the case of $H_i:\ \beta_i = \beta_j$, $1/A_i$ and $1/B_i$ are equal to zero. +- Equality constraints can be expressed as $$|\beta_i - \beta_j|\leq \xi$$ for small value of $\xi$. In this case, however, for vague encompassing prior $1/A_i$ is still virtually independent but $1/B_i$ is not. The more diffuse the encompassing prior is the smaller the Bayes factor in favor of the unconstrained model (similar to *Lindley’s paradox* where results are influenced by the vagueness of the prior distribution of the parameters involved). Sensitivity analysis is required considering both the choice of prior distributions and the value of $\xi$. + +#### Consideration {-} + +This is a possible approach. However, in our hypotheses there are several equality constraints. Thus, it would be important to properly evaluate them. + +### Inequality and Equality Constraints: A Second Approach {#in-eq-const} + +@guApproximatedAdjustedFractional2018 and @mulderBayesFactorTesting2019 presented an approach to evaluate **inequality and equality constraints** still based on the principle of encompassing priors. + +Informative hypotheses using equality and/or inequality constraints can be expressed as + +$$ +H_i:\ R_E\theta = r_E,\ R_I\theta > r_I +$$ +where $R_E$ and $R_I$ are the restriction matrices for equality and inequality constraints, respectively, and $r_E$ and $r_I$ contain constants. + +Following encompassing prior approach, given an unconstrained model $M_u$ and a nested model $M_i$, the Bayes Factor is computed as + +\begin{aligned} +BF_{iu} = \frac{Pr(\text{Inequality Const} | \text{Equality Const}, \text{Data})}{\pi_u(\text{Inequality Const} | \text{Equality Const})} &\frac{Pr(\text{Equality Const}| \text{Data})}{\pi_u( \text{Equality Const})} = \\ \\ \frac{Pr(R_I\theta > r_I | R_E\theta = r_E, Y)}{\pi_u(R_I\theta > r_I | R_E\theta = r_E)} &\frac{Pr(R_E\theta = r_E| Y)}{\pi_u( R_E\theta = r_E)}. +\end{aligned} + +Where, the first part is a ratio between conditional posterior and conditional prior **probabilities** that the inequality constraints hold under $M_u$. The second part is a ratio between marginal posterior and marginal prior **densities** of the equality constraints under $M_u$, the well-known *Savage–Dickey density ratio*. + +Note that the numerators provide information about the *fit* of $M_i$ relative to $M_u$ and the denominators provide information about the *complexity* of $M_i$ relative to $M_u$. + +#### Pro {-} + +- It allows to evaluate inequality and equality constraints. + +#### Cons {-} + +- The encompassing prior and posterior distributions are required in analytic form to compute conditional probabilities and marginal densities. + +#### Consideration {-} + +Obtaining the analytic posterior distribution is often difficult. To overcome this issue, @guApproximatedAdjustedFractional2018 and @mulderBayesFactorTesting2019 proposed two different solutions both based on **Normal approximations**. Let's consider the two approaches separately. + +#### @mulderBayesFactorTesting2019 {#mulder} + +Authors used to different solutions for computing the posterior and the prior: + +- **Computing Posterior** - a Metropolis-Hastings algorithm is used to obtain draws from the posterior of the unconstrained model. Subsequently the posterior was assumed to be approximately normally distributed +$$ +Pr(\theta|Y) \sim \mathcal{N}(\mu_\theta, \Sigma_\theta) +$$ +where $\mu_\theta$ and $\Sigma_\theta$ are the posterior mean and covariance matrix estimated from the posterior draws. Thus, conditional posterior probability and marginal posterior density are computed based on the normal approximation. + +- **Computing Prior** - in their case priors are not normally distributed. To estimate the marginal prior density they considered the fact that for a sufficiently small $\delta$, +$$ +Pr(|\theta_{E_1} - r_{E_1}| < \frac{\delta}{2}, \ldots, |\theta_{E_i} - r_{E_i}| < \frac{\delta}{2}|M_u) \approx \delta^i \times \pi_u(\theta_E = r_E) +$$ +Thus, the probability can be estimated considering the proportion of draws satisfying the constraints. To compute conditional prior probability, they used a normal approximation because this probability is not sensitive to the exact distributional form as long as they are symmetrical distributed around the value of the constraints, usually zero. The mean and covariance matrix were estimated based on the prior draws that satisfied the equality constraints. To do that they used the same approximation as before (i.e., $|\theta_{E_i} - r_{E_i}| < \frac{\delta}{2}$). + + + + +#### @guApproximatedAdjustedFractional2018 + +In this paper, the Approximated Adjusted Fractional Bayes Factors (AAFBF) method is presented. To avoid ad hoc or subjective specification of the unconstrained prior, the fractional Bayes factor approach is used. That is, proper default prior were defined updating non informative priors according to a a fraction *b* of the data +$$ +\text{Proper Prior: }\pi_u(\theta|Y_b) \approx Pr(Y_b|\theta) \times \pi_u(\theta). +$$ +The choice of the fraction $b$ also plays a crucial in the AAFBF and is extensively discussed in the article. + +Prior and posterior distributions are assumed to be approximately normally distributed due to the large-sample theory: +\begin{aligned} +\text{Prior: }& \pi_u(\theta|Y_b) \sim \mathcal{N}(\hat{\theta}, \hat{\Sigma}_\theta/b),\\ \\ +\text{Posterior: }& \pi_u(\theta|Y) \sim \mathcal{N}(\hat{\theta}, \hat{\Sigma}_\theta). +\end{aligned} +Where $\hat{\theta}$ and $\hat{\Sigma}_\theta$ denote the maximum likelihood estimate and covariance matrix of $\theta$, respectively. + +Moreover, authors proposed an adjustment of the prior mean for centering the prior distribution of $\theta$ on the boundary of the constrained space. This ensure that $Pr(\theta>0)$ and $Pr(\theta<0)$ are equal for symmetric distributions and follow common assumption that small effects are more likely a prior than large effects. To facilitate the computation of the adjusted prior mean, the parameters are transformed, +$$ +\beta = \left[\begin{array}{cc}\beta_E\\ \beta_I\end{array}\right] = \left[\begin{array}{cc}R_E \theta - r_E\\ R_I \theta - r_I\end{array}\right]. +$$ +In this way, informative hypothesis becomes $H_i:\ \beta_E=0, \beta_I>0$ and the prior mean vector is simply zero for the new parameter vector $\beta$. This parameter transformation from $\theta$ to $\beta$ simplifies the computation of the AAFBF. + +Note authors did not adjusted posterior mean since, according to large-sample theory, the prior has a negligible effect on the posterior for large samples. Thus, resulting prior and posterior distributions are: +\begin{aligned} +\text{Adjusted Prior: }& \pi_u^*(\beta|Y_b) \sim \mathcal{N}(0, \hat{\Sigma}_\beta/b),\\ \\ +\text{Posterior: }& \pi_u(\beta|Y) \sim \mathcal{N}(\hat{\beta}, \hat{\Sigma}_\beta). +\end{aligned} +Where $\hat{\beta} = R\hat{\theta}-r$ and $\hat{\Sigma}_\beta = R\hat{\Sigma}_\theta R^T$. + +### Present Study + +In our case, the unconstrained model is a Zero Inflated Poisson model with parameter vectors $\beta_p$ and $\beta_{\lambda}$ (see Section \@ref(model)): + +\begin{aligned} +y_i \sim& ZIPoisson(p_i, \lambda_i),\\ \\ +p_i =& \text{logit}^{-1}(X_i^T\beta_{p}),\\ \\ +\lambda_i =& \text{exp}(X_i^T\beta_{\lambda}). +\end{aligned} + +We can assume independent priors for each parameter $\beta_i \sim \mathcal{N}(0, \sigma)$, with $\sigma$ a large enough so the priors are vague. + +To evaluate informative hypothesis with equality and inequality constraints, we can follow the approach presented in Section \@ref(in-eq-const) considering the method proposed by @mulderBayesFactorTesting2019, see Section \@ref(mulder). + +In our case, both the prior and the posterior distributions are normally distributed (TODO: is it true?): +\begin{aligned} +\text{Prior: }& \Pr(\beta) \sim \mathcal{N}(0, \Sigma_{\sigma}),\\ \\ +\text{Posterior: }& \Pr(\beta|Y) \sim \mathcal{N}(\hat{\beta}, \hat{\Sigma}_{\beta}). +\end{aligned} +Where $\Sigma_{\sigma}$ is a diagonal matrix with all the diagonal elements equal to $\sigma$. Whereas, $\mu_\theta$ and $\Sigma_\theta$ are the posterior mean and covariance matrix estimated from the posterior draws. Posterior sample can be easily obtained in STAN. + +### Issues + +1. Ha senso questo approccio? +1. Come posso marginalizzare i valori delle densità? Vedi le loro formule (p.27 in basso). +1. Come ottenere le distribuzioni condizionate? Vedi le loro formule negli articoli (p.28). +1. Dovrei preferire invece l'Approximated Adjusted Fractional Bayes Factors method? +1. Nel modello ci saranno anche interazioni tra fattori, tutti i possibili gruppi. Devo pensare bene su come parametrizzare i contrasti. +1. Non sappiamo ancora quali saranno tutti i constraints e confronti ma sono tutti confronti tra gruppi. Devo essere sicuro che vengano rispettate le assunzioni delle ipotesi se uso l'Approximated Adjusted Fractional Bayes Factors method. Vedi "comparable hypotheses" formula 24 in @guApproximatedAdjustedFractional2018. +1. Potrei fare dell'hard coding codificando ogni gruppo con dummy variables. Non avrei "effetti principali" ma mi permette di evitare casini nel definire i confronti. (Nella codifica tradizionale i valori medi dei gruppi sarebbero dati da somme di più parametri, gli stessi parametri sarebbero usati per più gruppi e questo potrebbee compicare la definizione dei constraints sui parametri). +1. Devo stare attento che le prior non favoriscano alcun constraints e questo è problematico con gli equality constraints. (Al massimo faccio un po' di sensitivity sulle prior) + +## Estimating the Constrained Model + +Hypothesis testing will idicate us the most likely model according to the data. To evaluate the model parameters, howevere, we need to estimate the parameter posteriors. + +Parameter posteriors of the unconstrained model can be simply obtained using MCMC sampling implemented in STAN. However, obtaining parameter posteriors of models including inequality constraints is not straightforward. Here we describe two possible direction: the first one based on STAN and the second one based Gibbs sampling. + +### STAN Approach + +Considering the thread in the STAN forum *"How to give prior to ordered parameters, and also to simplex?"* ([link](https://discourse.mc-stan.org/t/how-to-give-prior-to-ordered-parameters-and-also-to-simplex/12041)), there is one answer ([link](https://discourse.mc-stan.org/t/how-to-give-prior-to-ordered-parameters-and-also-to-simplex/12041/7?u=claudio.zandonella)) suggesting a possible solution based on uniform priors. + +Given the restriction $a<\theta_1<\theta_2<\theta_3 \beta_j > \beta_k$ with $\beta_{i,j,k} \stackrel{i.i.d.}{\sim} \pi(\beta)$, the prior probability of the constraints is independent of the actual choice of the encompassing prior, and subsequently also $1/B_i$ is not influenced. + +#### Cons {-} + +- It is not possible to evaluate **equality constraints**. In the case of $H_i:\ \beta_i = \beta_j$, $1/A_i$ and $1/B_i$ are equal to zero. +- Equality constraints can be expressed as $$|\beta_i - \beta_j|\leq \xi$$ for small value of $\xi$. In this case, however, for vague encompassing prior $1/A_i$ is still virtually independent but $1/B_i$ is not. The more diffuse the encompassing prior is the smaller the Bayes factor in favor of the unconstrained model (similar to *Lindley’s paradox* where results are influenced by the vagueness of the prior distribution of the parameters involved). Sensitivity analysis is required considering both the choice of prior distributions and the value of $\xi$. + +#### Consideration {-} + +This is a possible approach. However, in our hypotheses there are several equality constraints. Thus, it would be important to properly evaluate them. + +## Inequality and Equality Constraints: A Second Approach {#in-eq-const} + +@guApproximatedAdjustedFractional2018 and @mulderBayesFactorTesting2019 presented an approach to evaluate **inequality and equality constraints** still based on the principle of encompassing priors. + +Informative hypotheses using equality and/or inequality constraints can be expressed as + +$$ +H_i:\ R_E\theta = r_E,\ R_I\theta > r_I +$$ +where $R_E$ and $R_I$ are the restriction matrices for equality and inequality constraints, respectively, and $r_E$ and $r_I$ contain constants. + +Following encompassing prior approach, given an unconstrained model $M_u$ and a nested model $M_i$, the Bayes Factor is computed as + +\begin{aligned} +BF_{iu} = \frac{Pr(\text{Inequality Const} | \text{Equality Const}, \text{Data})}{\pi_u(\text{Inequality Const} | \text{Equality Const})} &\frac{Pr(\text{Equality Const}| \text{Data})}{\pi_u( \text{Equality Const})} = \\ \\ \frac{Pr(R_I\theta > r_I | R_E\theta = r_E, Y)}{\pi_u(R_I\theta > r_I | R_E\theta = r_E)} &\frac{Pr(R_E\theta = r_E| Y)}{\pi_u( R_E\theta = r_E)}. +\end{aligned} + +Where, the first part is a ratio between conditional posterior and conditional prior **probabilities** that the inequality constraints hold under $M_u$. The second part is a ratio between marginal posterior and marginal prior **densities** of the equality constraints under $M_u$, the well-known *Savage–Dickey density ratio*. + +Note that the numerators provide information about the *fit* of $M_i$ relative to $M_u$ and the denominators provide information about the *complexity* of $M_i$ relative to $M_u$. + +#### Pro {-} + +- It allows to evaluate inequality and equality constraints. + +#### Cons {-} + +- The encompassing prior and posterior distributions are required in analytic form to compute conditional probabilities and marginal densities. + +#### Consideration {-} + +Obtaining the analytic posterior distribution is often difficult. To overcome this issue, @guApproximatedAdjustedFractional2018 and @mulderBayesFactorTesting2019 proposed two different solutions both based on **Normal approximations**. Let's consider the two approaches separately. + +### @mulderBayesFactorTesting2019 {#mulder} + +Authors used to different solutions for computing the posterior and the prior: + +- **Computing Posterior** - a Metropolis-Hastings algorithm is used to obtain draws from the posterior of the unconstrained model. Subsequently the posterior was assumed to be approximately normally distributed +$$ +Pr(\theta|Y) \sim \mathcal{N}(\mu_\theta, \Sigma_\theta) +$$ +where $\mu_\theta$ and $\Sigma_\theta$ are the posterior mean and covariance matrix estimated from the posterior draws. Thus, conditional posterior probability and marginal posterior density are computed based on the normal approximation. + +- **Computing Prior** - in their case priors are not normally distributed. To estimate the marginal prior density they considered the fact that for a sufficiently small $\delta$, +$$ +Pr(|\theta_{E_1} - r_{E_1}| < \frac{\delta}{2}, \ldots, |\theta_{E_i} - r_{E_i}| < \frac{\delta}{2}|M_u) \approx \delta^i \times \pi_u(\theta_E = r_E) +$$ +Thus, the probability can be estimated considering the proportion of draws satisfying the constraints. To compute conditional prior probability, they used a normal approximation because this probability is not sensitive to the exact distributional form as long as they are symmetrical distributed around the value of the constraints, usually zero. The mean and covariance matrix were estimated based on the prior draws that satisfied the equality constraints. To do that they used the same approximation as before (i.e., $|\theta_{E_i} - r_{E_i}| < \frac{\delta}{2}$). + + + + +### @guApproximatedAdjustedFractional2018 + +In this paper, the Approximated Adjusted Fractional Bayes Factors (AAFBF) method is presented. To avoid ad hoc or subjective specification of the unconstrained prior, the fractional Bayes factor approach is used. That is, proper default prior were defined updating non informative priors according to a a fraction *b* of the data +$$ +\text{Proper Prior: }\pi_u(\theta|Y_b) \approx Pr(Y_b|\theta) \times \pi_u(\theta). +$$ +The choice of the fraction $b$ also plays a crucial in the AAFBF and is extensively discussed in the article. + +Prior and posterior distributions are assumed to be approximately normally distributed due to the large-sample theory: +\begin{aligned} +\text{Prior: }& \pi_u(\theta|Y_b) \sim \mathcal{N}(\hat{\theta}, \hat{\Sigma}_\theta/b),\\ \\ +\text{Posterior: }& \pi_u(\theta|Y) \sim \mathcal{N}(\hat{\theta}, \hat{\Sigma}_\theta). +\end{aligned} +Where $\hat{\theta}$ and $\hat{\Sigma}_\theta$ denote the maximum likelihood estimate and covariance matrix of $\theta$, respectively. + +Moreover, authors proposed an adjustment of the prior mean for centering the prior distribution of $\theta$ on the boundary of the constrained space. This ensure that $Pr(\theta>0)$ and $Pr(\theta<0)$ are equal for symmetric distributions and follow common assumption that small effects are more likely a prior than large effects. To facilitate the computation of the adjusted prior mean, the parameters are transformed, +$$ +\beta = \left[\begin{array}{cc}\beta_E\\ \beta_I\end{array}\right] = \left[\begin{array}{cc}R_E \theta - r_E\\ R_I \theta - r_I\end{array}\right]. +$$ +In this way, informative hypothesis becomes $H_i:\ \beta_E=0, \beta_I>0$ and the prior mean vector is simply zero for the new parameter vector $\beta$. This parameter transformation from $\theta$ to $\beta$ simplifies the computation of the AAFBF. + +Note authors did not adjusted posterior mean since, according to large-sample theory, the prior has a negligible effect on the posterior for large samples. Thus, resulting prior and posterior distributions are: +\begin{aligned} +\text{Adjusted Prior: }& \pi_u^*(\beta|Y_b) \sim \mathcal{N}(0, \hat{\Sigma}_\beta/b),\\ \\ +\text{Posterior: }& \pi_u(\beta|Y) \sim \mathcal{N}(\hat{\beta}, \hat{\Sigma}_\beta). +\end{aligned} +Where $\hat{\beta} = R\hat{\theta}-r$ and $\hat{\Sigma}_\beta = R\hat{\Sigma}_\theta R^T$. + +## Present Study + +In our case, the unconstrained model is a Zero Inflated Poisson model with parameter vectors $\beta_p$ and $\beta_{\lambda}$ (see Section \@ref(model)): + +\begin{aligned} +y_i \sim& ZIPoisson(p_i, \lambda_i),\\ \\ +p_i =& \text{logit}^{-1}(X_i^T\beta_{p}),\\ \\ +\lambda_i =& \text{exp}(X_i^T\beta_{\lambda}). +\end{aligned} + +We can assume independent priors for each parameter $\beta_i \sim \mathcal{N}(0, \sigma)$, with $\sigma$ a large enough so the priors are vague. + +To evaluate informative hypothesis with equality and inequality constraints, we can follow the approach presented in Section \@ref(in-eq-const) considering the method proposed by @mulderBayesFactorTesting2019, see Section \@ref(mulder). + +In our case, both the prior and the posterior distributions are normally distributed (TODO: is it true?): +\begin{aligned} +\text{Prior: }& \Pr(\beta) \sim \mathcal{N}(0, \Sigma_{\sigma}),\\ \\ +\text{Posterior: }& \Pr(\beta|Y) \sim \mathcal{N}(\hat{\beta}, \hat{\Sigma}_{\beta}). +\end{aligned} +Where $\Sigma_{\sigma}$ is a diagonal matrix with all the diagonal elements equal to $\sigma$. Whereas, $\mu_\theta$ and $\Sigma_\theta$ are the posterior mean and covariance matrix estimated from the posterior draws. Posterior sample can be easily obtained in STAN. + +## Issues + +1. Ha senso questo approccio? +1. Come posso marginalizzare i valori delle densità? Vedi le loro formule (p.27 in basso). +1. Come ottenere le distribuzioni condizionate? Vedi le loro formule negli articoli (p.28). +1. Dovrei preferire invece l'Approximated Adjusted Fractional Bayes Factors method? +1. Nel modello ci saranno anche interazioni tra fattori, tutti i possibili gruppi. Devo pensare bene su come parametrizzare i contrasti. +1. Non sappiamo ancora quali saranno tutti i constraints e confronti ma sono tutti confronti tra gruppi. Devo essere sicuro che vengano rispettate le assunzioni delle ipotesi se uso l'Approximated Adjusted Fractional Bayes Factors method. Vedi "comparable hypotheses" formula 24 in @guApproximatedAdjustedFractional2018. +1. Potrei fare dell'hard coding codificando ogni gruppo con dummy variables. Non avrei "effetti principali" ma mi permette di evitare casini nel definire i confronti. (Nella codifica tradizionale i valori medi dei gruppi sarebbero dati da somme di più parametri, gli stessi parametri sarebbero usati per più gruppi e questo potrebbee compicare la definizione dei constraints sui parametri). +1. Devo stare attento che le prior non favoriscano alcun constraints e questo è problematico con gli equality constraints. (Al massimo faccio un po' di sensitivity sulle prior) + +# Estimating the Constrained Model + +Hypothesis testing will idicate us the most likely model according to the data. To evaluate the model parameters, howevere, we need to estimate the parameter posteriors. + +Parameter posteriors of the unconstrained model can be simply obtained using MCMC sampling implemented in STAN. However, obtaining parameter posteriors of models including inequality constraints is not straightforward. Here we describe two possible direction: the first one based on STAN and the second one based Gibbs sampling. + +## STAN Approach + +Considering the thread in the STAN forum *"How to give prior to ordered parameters, and also to simplex?"* ([link](https://discourse.mc-stan.org/t/how-to-give-prior-to-ordered-parameters-and-also-to-simplex/12041)), there is one answer ([link](https://discourse.mc-stan.org/t/how-to-give-prior-to-ordered-parameters-and-also-to-simplex/12041/7?u=claudio.zandonella)) suggesting a possible solution based on uniform priors. + +Given the restriction $a<\theta_1<\theta_2<\theta_3 + diff --git a/Documents/Report-Project/bibliography.bib b/Documents/Report-Project/bibliography.bib new file mode 100644 index 0000000..108da06 --- /dev/null +++ b/Documents/Report-Project/bibliography.bib @@ -0,0 +1,18 @@ + +@article{fisherFrequencyDistributionValues1915, + title = {Frequency {{Distribution}} of the {{Values}} of the {{Correlation Coefficient}} in {{Samples}} from an {{Indefinitely Large Population}}}, + author = {Fisher, R. A.}, + date = {1915}, + journaltitle = {Biometrika}, + volume = {10}, + pages = {507}, + issn = {00063444}, + doi = {10.2307/2331838}, + eprint = {2331838}, + eprinttype = {jstor}, + file = {/Users/claudio/MEGA/Zotero/Fisher_1915_Frequency Distribution of the Values of the Correlation Coefficient in Samples.pdf}, + langid = {english}, + number = {4} +} + + diff --git a/Documents/Report-Project/css/custom.css b/Documents/Report-Project/css/custom.css new file mode 100644 index 0000000..0c43243 --- /dev/null +++ b/Documents/Report-Project/css/custom.css @@ -0,0 +1,27 @@ +/*----------LOGO above TOC---------*/ + +#TOC::before { + content: ""; + display: block; + height: 200px; + margin: 2em 10px 20px 10px; + background-image: url("hex_psicostat.png"); + background-size: contain; + background-position: center center; + background-repeat: no-repeat; +} + +/*----------Color active TOC---------*/ +.list-group-item.active, .list-group-item.active:focus, .list-group-item.active:hover { + background-color: #8C3031; + border-color: #8C3031; +} + +/*----------Color links---------*/ +a { + color: #8C3031; + font-weight: bold; +} + +/*----------line spacing---------*/ +p {line-height: 1.5em;} diff --git a/Documents/Report-Project/css/hex_psicostat.png b/Documents/Report-Project/css/hex_psicostat.png new file mode 100644 index 0000000..1b36060 Binary files /dev/null and b/Documents/Report-Project/css/hex_psicostat.png differ diff --git a/Documents/Report-Project/figures/hex_psicostat.png b/Documents/Report-Project/figures/hex_psicostat.png new file mode 100644 index 0000000..1b36060 Binary files /dev/null and b/Documents/Report-Project/figures/hex_psicostat.png differ diff --git a/Documents/Report-Project/figures/logo_psicostat.png b/Documents/Report-Project/figures/logo_psicostat.png new file mode 100644 index 0000000..874a682 Binary files /dev/null and b/Documents/Report-Project/figures/logo_psicostat.png differ diff --git a/Documents/analisi_d1.Rmd b/Documents/analisi_d1.Rmd new file mode 100644 index 0000000..12b8064 --- /dev/null +++ b/Documents/analisi_d1.Rmd @@ -0,0 +1,376 @@ +--- +title: "Analisi dei Cluster - Zandonella Callegher et al" +author: '' +date: "Maggio 2021" +output: + word_document: + toc: yes + toc_depth: '3' + html_notebook: + number_sections: no + toc: yes + toc_depth: 3 + toc_float: yes + html_document: + number_sections: no + toc: yes + toc_depth: 3 + toc_float: yes +editor_options: + chunk_output_type: console +--- + + +```{r,warning=FALSE} + +rm(list=ls()) +options(digits = 3) +``` + +Librerie +```{r,warning=FALSE,include=TRUE} +library(car) +library(semPlot) +library(foreign) +library(lavaan) +library(psych) +library(effects) +library(sciplot) +library(psych) +``` + +```{r,include=FALSE} +# effect size +etas=function(m.lm){ + library(car) + ris=NULL + anovam=Anova(m.lm,type="III") + l.SS=length(anovam$"Sum Sq") + SStotal=sum(anovam$"Sum Sq"[2:l.SS]) + eta.q=(anovam$"Sum Sq"[2:(l.SS-1)])/SStotal + peta.q=(anovam$"Sum Sq"[2:(l.SS-1)]) / ( (anovam$"Sum Sq"[2:(l.SS-1)]) + (anovam$"Sum Sq"[l.SS])) + ris=data.frame(eta.q,peta.q) + row.names(ris)=row.names(anovam)[2:(l.SS-1)] + return(ris) +} +``` + +```{r,warning=FALSE,include=TRUE} +####################################################################################### +# # +setwd("~/Documents/Mega/Paper/In preparation/Zandonella et al/Analisi_data_d1") + +load("/Users/tatiana/Documents/Mega/Paper/In preparation/Zandonella et al/Analisi_data_d1/ddef_datiECR.rda") + +``` + + + +```{r} +ddef=ddef_datiECR[ddef_datiECR$age_year<12.30,] +ddef<- subset(ddef, ddef$fas !="NA") +nrow(ddef) #847 +save(ddef,file=paste('ddef.rda',sep="")) +``` + +```{r,warning=FALSE} +load("/Users/tatiana/Documents/Mega/Paper/In preparation/Zandonella et al/Analisi_data_d1/ddef.rda") +``` + + + +```{r,warning=FALSE} +####################################################################################### +######################### DISTRIBUZIONE SDQ TOTALE +# + +#quantile(ddef$ebdtot_somma, prob = c(0.00, 0.30, 0.40, 0.70, 1), na.rm=TRUE) + +#cut-off: +#0-13 normale +#borderline 14-16 +#problematico 17-40 +ddef$SDQgroup <-ifelse(ddef$sdqtot_sum < 14, "1", + ifelse((ddef$sdqtot_sum >= 14) & (ddef$sdqtot_sum<=16),"2","3")) + +prop.table(table(ddef$SDQgroup)) + +# pochi problemi! +``` + + + +```{r} +#se dovesse servire lavorare su dati primaria +primaria<-ddef[ddef$ordine_scolastico=="PRIMARIA",] +save(primaria,file=paste('primaria.rda',sep="")) +``` + + + +#Cluster mamma +```{r,warning=FALSE} +####################################################################################### +#1. ######################## CLUSTER +# +dcM <- dist(ddef[,c(60:71)], method = "euclidean") # distance matrix +fit <- hclust(dcM, method="ward.D") +plot(fit, main = "Dendrogramma") +# dal dendrogramma sembrano emergere 4 cluster +rect.hclust(fit, k=4, border="red") +ddef$gruppiattM <- as.factor(cutree(fit, k=4)) + +save(ddef,file=paste('ddef.rda',sep="")) + +####################################################### +table(ddef$gruppiattM) +prop.table(table(ddef$gruppiattM)) +prop.table(table(ddef$gruppiattM,ddef$genere),2) +summary(table(ddef$gruppiattM,ddef$genere)) +summary(table(ddef$gruppiattM,ddef$fas)) + +``` + +### Profilo grafico cluster +```{r} +fattore=rep(c("Ansia","Evitamento"),rep(dim(ddef)[1],2)) +y=c(scale(ddef$Anxm),scale(ddef$Avm)) +y1=c(ddef$Anxm,ddef$Avm) +gruppo=rep(ddef$gruppiattM,2) +dprof=data.frame(y,y1,fattore,gruppo) + + +dprof$nomi.gruppi=factor(dprof$gruppo,labels=c("Gruppo1: Sicuri (n = 231)","Gruppo 2: Ansiosi (n = 286)","Gruppo 3: Ansiosi-Evitanti (n = 100)","Gruppo 4: Evitanti (n = 230)")) + + + +``` + +```{r} +par(mar=c(6, 5, 4, 2) + 0.1) + +bargraph.CI(fattore, y, group = nomi.gruppi, data=dprof,legend=T,ylim=c(-2.4,1.6),x.leg=1.2, y.leg=-1.2,ylab="Punteggi standardizzati",xlab="",cex.lab=.8,cex.leg=.8, col=c("yellow","green","blue","red")) + +save(ddef,file=paste('ddef.rda',sep="")) +``` + + + +### Valori medi per i 4 gruppi +```{r,warning=FALSE} +tapply(ddef$Anxm,ddef$gruppiattM,mean) +tapply(ddef$Avm,ddef$gruppiattM,mean) + + +tapply(ddef$Anxm,ddef$gruppiattM,sd) +tapply(ddef$Avm,ddef$gruppiattM,sd) +``` + + +### Inferenza sui 4 gruppi rispetto a attaccamento mamma +```{r,warning=FALSE} +# mamma ansia +mod.anx=lm(Anxm~gruppiattM,ddef) +summary(mod.anx) +Anova(mod.anx) +etas(mod.anx) + + +# mamma evitamento +mod.av=lm(Avm~gruppiattM,ddef) +summary(mod.av) +Anova(mod.av) +etas(mod.av) + +``` + + +### Prove su internalizing and externalizing + +```{r} +m1.1<-lm(externalizing~gruppiattM,ddef) +Anova(m1.1) +summary(m1.1) + + +g1<- effect("gruppiattM ", m1.1) +plot(g1,ylab="problemi", xlab="gruppi attaccamento mamma",main="",ci.stile="band",band.colors="red") + +# R-squared basso. A livello descrittivo risultati interessanti e corrispondenti alle attese. Gruppo 3 (alta ansia, alto evitamento) maggiormente problematico + + +m1.2<-lm(externalizing~gruppiattM*genere+fas,ddef) +summary(m1.2) +Anova(m1.2) + +g1.2<- effect("gruppiattM:genere ", m1.2) +plot(g1.2,ylab="problemi", xlab="gruppi attaccamento mamma",main="",ci.stile="band",band.colors="red") + +#no interazione con genere +``` + +```{r} +m2.1<-lm(internalizing~gruppiattM,ddef) +summary(m2.1) +Anova(m2.1) + + + +g2.1<- effect("gruppiattM ", m2.1) +plot(g2.1,ylab="problemi", xlab="gruppi attaccamento mamma",main="",ci.stile="band",band.colors="red") + +# R-squared basso. A livello descrittivo risultati interessanti e corrispondenti alle attese. Gruppo 3 (alta ansia, alto evitamento) maggiormente problematico + + +m2.2<-lm(internalizing~gruppiattM*genere,ddef) +summary(m2.2) +Anova(m2.2) + + + +g2.2<- effect("gruppiattM:genere ", m2.2) +plot(g2.2,ylab="problemi", xlab="gruppi attaccamento mamma",main="",ci.stile="band",band.colors="red") + +# no effetto genere ma a livello descrittivo risultati interessanti. vedere in particolare gruppo 3 femmine + +``` + + +#Cluster papà +```{r,warning=FALSE} +####################################################################################### +#1. ######################## CLUSTER +# +dcP <- dist(ddef[,c(99:110)], method = "euclidean") # distance matrix +fit <- hclust(dcP, method="ward.D") +plot(fit, main = "Dendrogramma") +# dal dendrogramma sembrano emergere 4 cluster +rect.hclust(fit, k=4, border="red") +ddef$gruppiattP <- as.factor(cutree(fit, k=4)) + +save(ddef,file=paste('ddef.rda',sep="")) + + +table(ddef$gruppiattP) +prop.table(table(ddef$gruppiattP)) +prop.table(table(ddef$gruppiattP,ddef$genere),2) +summary(table(ddef$gruppiattP,ddef$genere)) +summary(table(ddef$gruppiattP,ddef$fas)) + +``` + +### Profilo grafico cluster +```{r} +fattore=rep(c("Ansia","Evitamento"),rep(dim(ddef)[1],2)) +y=c(scale(ddef$Anxp),scale(ddef$Avp)) +y1=c(ddef$Anxp,ddef$Avp) +gruppo=rep(ddef$gruppiattP,2) +dprof=data.frame(y,y1,fattore,gruppo) +``` + + +### Valori medi per i 4 gruppi +```{r,warning=FALSE} +tapply(ddef$Anxp,ddef$gruppiattP,mean) +tapply(ddef$Avp,ddef$gruppiattP,mean) + + +tapply(ddef$Anxp,ddef$gruppiattP,sd) +tapply(ddef$Avp,ddef$gruppiattP,sd) +``` + +### Inferenza sui 4 gruppi rispetto a attaccamento papà +```{r,warning=FALSE} +# papà ansia +mod.anx=lm(Anxp~gruppiattP,ddef) +summary(mod.anx) +Anova(mod.anx) # +etas(mod.anx) + + +# papà evitamento +mod.av=lm(Avp~gruppiattP,ddef) +summary(mod.av) +Anova(mod.av) +etas(mod.av) + +``` + + + +```{r,warning=FALSE} + +dprof$nomi.gruppi=factor(dprof$gruppo,labels=c("Gruppo1: Sicuri (n = 206)","Gruppo 2: Ansiosi (n = 230)","Gruppo 3: Evitanti (n = 311 )","Gruppo 4: Ansiosi-Evitanti (n = 100)")) + + +par(mar=c(6, 5, 4, 2) + 0.1) + +bargraph.CI(fattore, y, group = nomi.gruppi, data=dprof,legend=T,ylim=c(-2.4,1.6),x.leg=1.2, y.leg=-1.2,ylab="Punteggi standardizzati",xlab="",cex.lab=.8,cex.leg=.8, col=c("yellow","green","blue","red")) + +#N.B: nel gruppo 2 anche se nominati ansiosi i livelli sono comunque bassi +box() + +save(ddef,file=paste('ddef.rda',sep="")) +``` + + +###Prove su internalizing e externalizing +```{r} +m3.1<-lm(externalizing~gruppiattP,ddef) +Anova(m3.1) +summary(m3.1) + + +g3.1<- effect("gruppiattP ", m3.1) +plot(g3.1,ylab="problemi esternalizzanti", xlab="gruppi attaccamento papà",main="",ci.stile="band",band.colors="red") + +# R-squared bassissimo. A livello descrittivo Gruppo 4 (alta ansia, alto evitamento) maggiormente problematico (come gruppo 3 mamma corrispondente a alta ansia alto evitamento) + + +m3.2<-lm(externalizing~gruppiattP*genere+fas,ddef) +summary(m3.2) +Anova(m3.2) + +g3.2<- effect("gruppiattP:genere ", m3.2) +plot(g3.2,ylab="problemi esternalizzanti", xlab="gruppi attaccamento papà",main="",ci.stile="band",band.colors="red") + +#no interazione con genere. In generale maschi più problemi esternalizzanti (risultati coerenti con la letteratura) +``` + +```{r} +m4.1<-lm(internalizing~gruppiattP,ddef) +summary(m4.1) +Anova(m4.1) + + + +g4.1<- effect("gruppiattP ", m4.1) +plot(g4.1,ylab="problemi", xlab="gruppi attaccamento papà",main="",ci.stile="band",band.colors="red") + +# R-squared bassissimo. A livello descrittivo risultati interessanti e corrispondenti alle attese. Gruppo 4 (alta ansia, alto evitamento) maggiormente problematico + + +m4.2<-lm(internalizing~gruppiattP*genere,ddef) +summary(m4.2) +Anova(m4.2) + + + +g4.2<- effect("gruppiattP:genere ", m4.2) +plot(g4.2,ylab="problemi", xlab="gruppi attaccamento papà",main="",ci.stile="band",band.colors="red") + +# no effetto genere ma a livello descrittivo risultati interessanti. + +``` + +Note label gruppi: + +Gruoo1 mamma = Gruppo 1 papà = sicuri + +Gruoo2 mamma = Gruppo 2 papà = ansiosi + +Gruoo3 mamma = Gruppo 4 papà = ansiosi-evitanti + +Gruoo4 mamma = Gruppo 3 papà = evitanti + + diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..6ae9268 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,2 @@ +# Generated by roxygen2: do not edit by hand + diff --git a/R/Clusters.R b/R/Clusters.R new file mode 100644 index 0000000..129277a --- /dev/null +++ b/R/Clusters.R @@ -0,0 +1,76 @@ +#================================# +#==== Cluster Analysis ====# +#================================# + +#----- cluster_tatiana ---- + +cluster_tatiana <- function(data, parent = c("mother", "father")){ + parent <- match.arg(parent) + + if(parent == "mother"){ + cols_num <- 60:71 + } else { + cols_num <- 99:110 + } + dist_matrix <- dist(data[,cols_num], method = "euclidean") + fit_cluster <- hclust(dist_matrix, method="ward.D") + + return(fit_cluster) +} + +#---- get_cluster_var ---- + +get_cluster_var <- function(fit_cluster){ + as.factor(cutree(fit_cluster, k=4)) +} + +#----- plot_scores ---- + +plot_scores_mother <- function(data){ + fattore=rep(c("Ansia","Evitamento"),rep(dim(data)[1],2)) + y=c(scale(data$Anxm),scale(data$Avm)) + y1=c(data$Anxm,data$Avm) + gruppo=rep(data$cluster_mother_tat,2) + dprof=data.frame(y,y1,fattore,gruppo) + + + dprof$nomi.gruppi=factor(dprof$gruppo,labels=c("Gruppo1: Sicuri (n = 231)","Gruppo 2: Ansiosi (n = 286)","Gruppo 3: Ansiosi-Evitanti (n = 100)","Gruppo 4: Evitanti (n = 230)")) + sciplot::bargraph.CI(fattore, y, group = nomi.gruppi, data=dprof,legend=T,ylim=c(-2.4,1.6),x.leg=1.2, y.leg=-1.2,ylab="Punteggi standardizzati",xlab="",cex.lab=.8,cex.leg=.8, col=c("yellow","green","blue","red")) +} + +plot_scores_father <- function(data){ + fattore=rep(c("Ansia","Evitamento"),rep(dim(data)[1],2)) + y=c(scale(data$Anxp),scale(data$Avp)) + y1=c(data$Anxp,data$Avp) + gruppo=rep(data$cluster_father_tat,2) + dprof=data.frame(y,y1,fattore,gruppo) + + dprof$nomi.gruppi=factor(dprof$gruppo,labels=c("Gruppo1: Sicuri (n = 206)","Gruppo 2: Ansiosi (n = 230)","Gruppo 3: Evitanti (n = 311 )","Gruppo 4: Ansiosi-Evitanti (n = 100)")) + sciplot::bargraph.CI(fattore, y, group = nomi.gruppi, data=dprof,legend=T,ylim=c(-2.4,1.6),x.leg=1.2, y.leg=-1.2,ylab="Punteggi standardizzati",xlab="",cex.lab=.8,cex.leg=.8, col=c("yellow","green","blue","red")) +} + +#---- mclust_BIC ---- + +mclust_BIC <- function(data, parent = c("mother", "father")){ + parent <- match.arg(parent) + + if(parent == "mother"){ + cols_num <- 60:71 + class <- data$cluster_mother_tat + } else { + cols_num <- 99:110 + class <- data$cluster_father_tat + } + + data_no_NA <- data[, cols_num] + data_no_NA[is.na(data_no_NA)] <- 1 + + X <- data_no_NA + + res <- mclust::mclustBIC(X) + + return(res) +} +#----- + + diff --git a/R/Models.R b/R/Models.R new file mode 100644 index 0000000..405592c --- /dev/null +++ b/R/Models.R @@ -0,0 +1,40 @@ +#================================# +#==== Model Definition ====# +#================================# + +#----- zip_brms ---- + +zip_brms <- function(data){ + + fit <- brms::brm(brms::bf( + internalizing_sum ~ genere + cluster_mother_tat + cluster_father_tat, + zi ~ genere + cluster_mother_tat + cluster_father_tat), + data = data, family = brms::zero_inflated_poisson()) + + return(fit) +} + +#---- make_stan_data ---- + +make_stan_data <- function(data = data_cluster_tat){ + formula <- brms::bf( + internalizing_sum ~ genere + cluster_mother_tat + cluster_father_tat, + zi ~ genere + cluster_mother_tat + cluster_father_tat) + + formula <- brms:::validate_formula(formula, data = data, family = brms::zero_inflated_poisson(), + autocor = NULL, sparse = NULL, cov_ranef = NULL) + + bterms <- brms:::brmsterms(formula) + + + data_name <- brms:::substitute_name(data) + data <- brms:::validate_data(data, bterms = bterms, data2 = NULL, knots = NULL) + attr(data, "data_name") <- data_name + + sdata <- brms:::.make_standata(bterms, data = data, prior = NULL, + data2 = NULL, stanvars = NULL, threads = NULL) + + return(sdata) +} + +#----- diff --git a/R/Munge_data.R b/R/Munge_data.R new file mode 100644 index 0000000..f9f4719 --- /dev/null +++ b/R/Munge_data.R @@ -0,0 +1,22 @@ +#============================# +#==== Data Munging ====# +#============================# + +#----- munge_data ---- + +munge_data <- function(data){ + data <- data[data$age_year<12.30,] + data <- subset(data, data$fas !="NA") + + data <- data %>% + mutate_at(c("externalizing_sum", "internalizing_sum"), round, 0) + + return(data) +} + + +#----- + + + + diff --git a/R/Plan.R b/R/Plan.R new file mode 100644 index 0000000..683cba6 --- /dev/null +++ b/R/Plan.R @@ -0,0 +1,43 @@ +############################# +#### Analysis Plan #### +############################# + + +#---- get_analysis_plan ---- + +#' @title Plan of the Analysis +#' +#' @description Plan of the Analysis +#' +#' @return a \code{drake_plan} object of class \code{tbl_df} with the plan of +#' the analysis +#' + +get_analysis_plan <- function(){ + drake::drake_plan( + + #---- Munge Data ---- + raw_data = readRDS("Data/data_ECR.rds"), + data_munged = munge_data(data = raw_data), + + + #---- Cluster Analysis ---- + cluster_mother_tat = cluster_tatiana(data = data_munged, parent = "mother"), + cluster_father_tat = cluster_tatiana(data = data_munged, parent = "father"), + + data_cluster_tat = cbind(data_munged, + cluster_mother_tat = get_cluster_var(cluster_mother_tat), + cluster_father_tat = get_cluster_var(cluster_father_tat)), + + # mclust check + mclust_mother = mclust_BIC(data = data_cluster_tat, parent = "mother"), + mclust_father = mclust_BIC(data = data_cluster_tat, parent = "father"), + + + #---- brms Models ---- + fit_int = zip_brms(data_cluster_tat), + stan_data = make_stan_data(data = data_cluster_tat) + ) +} + +#---- diff --git a/Stan/ZIP-model-prior.stan b/Stan/ZIP-model-prior.stan new file mode 100644 index 0000000..bc2b3f8 --- /dev/null +++ b/Stan/ZIP-model-prior.stan @@ -0,0 +1,153 @@ +// generated with brms 2.15.0 +functions { + /* zero-inflated poisson log-PDF of a single response + * Args: + * y: the response value + * lambda: mean parameter of the poisson distribution + * zi: zero-inflation probability + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_lpmf(int y, real lambda, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_lpmf(1 | zi), + bernoulli_lpmf(0 | zi) + + poisson_lpmf(0 | lambda)); + } else { + return bernoulli_lpmf(0 | zi) + + poisson_lpmf(y | lambda); + } + } + /* zero-inflated poisson log-PDF of a single response + * logit parameterization of the zero-inflation part + * Args: + * y: the response value + * lambda: mean parameter of the poisson distribution + * zi: linear predictor for zero-inflation part + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_logit_lpmf(int y, real lambda, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_logit_lpmf(1 | zi), + bernoulli_logit_lpmf(0 | zi) + + poisson_lpmf(0 | lambda)); + } else { + return bernoulli_logit_lpmf(0 | zi) + + poisson_lpmf(y | lambda); + } + } + /* zero-inflated poisson log-PDF of a single response + * log parameterization for the poisson part + * Args: + * y: the response value + * eta: linear predictor for poisson distribution + * zi: zero-inflation probability + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_log_lpmf(int y, real eta, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_lpmf(1 | zi), + bernoulli_lpmf(0 | zi) + + poisson_log_lpmf(0 | eta)); + } else { + return bernoulli_lpmf(0 | zi) + + poisson_log_lpmf(y | eta); + } + } + /* zero-inflated poisson log-PDF of a single response + * log parameterization for the poisson part + * logit parameterization of the zero-inflation part + * Args: + * y: the response value + * eta: linear predictor for poisson distribution + * zi: linear predictor for zero-inflation part + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_log_logit_lpmf(int y, real eta, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_logit_lpmf(1 | zi), + bernoulli_logit_lpmf(0 | zi) + + poisson_log_lpmf(0 | eta)); + } else { + return bernoulli_logit_lpmf(0 | zi) + + poisson_log_lpmf(y | eta); + } + } + // zero-inflated poisson log-CCDF and log-CDF functions + real zero_inflated_poisson_lccdf(int y, real lambda, real zi) { + return bernoulli_lpmf(0 | zi) + poisson_lccdf(y | lambda); + } + real zero_inflated_poisson_lcdf(int y, real lambda, real zi) { + return log1m_exp(zero_inflated_poisson_lccdf(y | lambda, zi)); + } +} +data { + int N; // total number of observations + int Y[N]; // response variable + int K; // number of population-level effects + matrix[N, K] X; // population-level design matrix + int K_zi; // number of population-level effects + matrix[N, K_zi] X_zi; // population-level design matrix + int prior_only; // should the likelihood be ignored? +} +transformed data { + int Kc = K - 1; + matrix[N, Kc] Xc; // centered version of X without an intercept + vector[Kc] means_X; // column means of X before centering + int Kc_zi = K_zi - 1; + matrix[N, Kc_zi] Xc_zi; // centered version of X_zi without an intercept + vector[Kc_zi] means_X_zi; // column means of X_zi before centering + for (i in 2:K) { + means_X[i - 1] = mean(X[, i]); + Xc[, i - 1] = X[, i] - means_X[i - 1]; + } + for (i in 2:K_zi) { + means_X_zi[i - 1] = mean(X_zi[, i]); + Xc_zi[, i - 1] = X_zi[, i] - means_X_zi[i - 1]; + } +} +parameters { + // lambda regression parameters + real b_gender; + vector[3] db_mother; // order constraints mother + vector[3] b_father; // parameters farther + real Intercept; // temporary intercept for centered predictors + + // zi regression parameters + vector[Kc_zi] b_zi; // population-level effects + real Intercept_zi; // temporary intercept for centered predictors +} +transformed parameters { + vector[Kc] b; // population-level effects + + b[1] = b_gender; + b[2] = db_mother[1]; + b[3] = b[2] + db_mother[2]; + b[4] = db_mother[3]; + b[5:7] = b_father; +} +model { + // likelihood including constants + if (!prior_only) { + // initialize linear predictor term + vector[N] mu = Intercept + Xc * b; + // initialize linear predictor term + vector[N] zi = Intercept_zi + Xc_zi * b_zi; + for (n in 1:N) { + target += zero_inflated_poisson_log_logit_lpmf(Y[n] | mu[n], zi[n]); + } + } + // priors including constants + target += uniform_lpdf(db_mother[2] | 0, 5 - b[2]); + target += student_t_lpdf(Intercept | 3, 0.7, 2.5); + target += logistic_lpdf(Intercept_zi | 0, 1); +} +generated quantities { + // actual population-level intercept + real b_Intercept = Intercept - dot_product(means_X, b); + // actual population-level intercept + real b_zi_Intercept = Intercept_zi - dot_product(means_X_zi, b_zi); +} diff --git a/Stan/ZIP-model.stan b/Stan/ZIP-model.stan new file mode 100644 index 0000000..2411127 --- /dev/null +++ b/Stan/ZIP-model.stan @@ -0,0 +1,140 @@ +// generated with brms 2.15.0 +functions { + /* zero-inflated poisson log-PDF of a single response + * Args: + * y: the response value + * lambda: mean parameter of the poisson distribution + * zi: zero-inflation probability + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_lpmf(int y, real lambda, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_lpmf(1 | zi), + bernoulli_lpmf(0 | zi) + + poisson_lpmf(0 | lambda)); + } else { + return bernoulli_lpmf(0 | zi) + + poisson_lpmf(y | lambda); + } + } + /* zero-inflated poisson log-PDF of a single response + * logit parameterization of the zero-inflation part + * Args: + * y: the response value + * lambda: mean parameter of the poisson distribution + * zi: linear predictor for zero-inflation part + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_logit_lpmf(int y, real lambda, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_logit_lpmf(1 | zi), + bernoulli_logit_lpmf(0 | zi) + + poisson_lpmf(0 | lambda)); + } else { + return bernoulli_logit_lpmf(0 | zi) + + poisson_lpmf(y | lambda); + } + } + /* zero-inflated poisson log-PDF of a single response + * log parameterization for the poisson part + * Args: + * y: the response value + * eta: linear predictor for poisson distribution + * zi: zero-inflation probability + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_log_lpmf(int y, real eta, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_lpmf(1 | zi), + bernoulli_lpmf(0 | zi) + + poisson_log_lpmf(0 | eta)); + } else { + return bernoulli_lpmf(0 | zi) + + poisson_log_lpmf(y | eta); + } + } + /* zero-inflated poisson log-PDF of a single response + * log parameterization for the poisson part + * logit parameterization of the zero-inflation part + * Args: + * y: the response value + * eta: linear predictor for poisson distribution + * zi: linear predictor for zero-inflation part + * Returns: + * a scalar to be added to the log posterior + */ + real zero_inflated_poisson_log_logit_lpmf(int y, real eta, real zi) { + if (y == 0) { + return log_sum_exp(bernoulli_logit_lpmf(1 | zi), + bernoulli_logit_lpmf(0 | zi) + + poisson_log_lpmf(0 | eta)); + } else { + return bernoulli_logit_lpmf(0 | zi) + + poisson_log_lpmf(y | eta); + } + } + // zero-inflated poisson log-CCDF and log-CDF functions + real zero_inflated_poisson_lccdf(int y, real lambda, real zi) { + return bernoulli_lpmf(0 | zi) + poisson_lccdf(y | lambda); + } + real zero_inflated_poisson_lcdf(int y, real lambda, real zi) { + return log1m_exp(zero_inflated_poisson_lccdf(y | lambda, zi)); + } +} +data { + int N; // total number of observations + int Y[N]; // response variable + int K; // number of population-level effects + matrix[N, K] X; // population-level design matrix + int K_zi; // number of population-level effects + matrix[N, K_zi] X_zi; // population-level design matrix + int prior_only; // should the likelihood be ignored? +} +transformed data { + int Kc = K - 1; + matrix[N, Kc] Xc; // centered version of X without an intercept + vector[Kc] means_X; // column means of X before centering + int Kc_zi = K_zi - 1; + matrix[N, Kc_zi] Xc_zi; // centered version of X_zi without an intercept + vector[Kc_zi] means_X_zi; // column means of X_zi before centering + for (i in 2:K) { + means_X[i - 1] = mean(X[, i]); + Xc[, i - 1] = X[, i] - means_X[i - 1]; + } + for (i in 2:K_zi) { + means_X_zi[i - 1] = mean(X_zi[, i]); + Xc_zi[, i - 1] = X_zi[, i] - means_X_zi[i - 1]; + } +} +parameters { + vector[Kc] b; // population-level effects + real Intercept; // temporary intercept for centered predictors + vector[Kc_zi] b_zi; // population-level effects + real Intercept_zi; // temporary intercept for centered predictors +} +transformed parameters { +} +model { + // likelihood including constants + if (!prior_only) { + // initialize linear predictor term + vector[N] mu = Intercept + Xc * b; + // initialize linear predictor term + vector[N] zi = Intercept_zi + Xc_zi * b_zi; + for (n in 1:N) { + target += zero_inflated_poisson_log_logit_lpmf(Y[n] | mu[n], zi[n]); + } + } + // priors including constants + target += student_t_lpdf(Intercept | 3, 0.7, 2.5); + target += logistic_lpdf(Intercept_zi | 0, 1); +} +generated quantities { + // actual population-level intercept + real b_Intercept = Intercept - dot_product(means_X, b); + // actual population-level intercept + real b_zi_Intercept = Intercept_zi - dot_product(means_X_zi, b_zi); +} diff --git a/renv.lock b/renv.lock new file mode 100644 index 0000000..fab7cd5 --- /dev/null +++ b/renv.lock @@ -0,0 +1,1175 @@ +{ + "R": { + "Version": "4.1.0", + "Repositories": [ + { + "Name": "CRAN", + "URL": "https://cran.rstudio.com" + } + ] + }, + "Packages": { + "BH": { + "Package": "BH", + "Version": "1.75.0-0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e4c04affc2cac20c8fec18385cd14691" + }, + "Brobdingnag": { + "Package": "Brobdingnag", + "Version": "1.2-6", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "365629a3ac7243df8ed6ac88a606d9fe" + }, + "DBI": { + "Package": "DBI", + "Version": "1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "030aaec5bc6553f35347cbb1e70b1a17" + }, + "DT": { + "Package": "DT", + "Version": "0.18", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a7d6660c869d4f41f856504828af4645" + }, + "MASS": { + "Package": "MASS", + "Version": "7.3-54", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0e59129db205112e3963904db67fd0dc" + }, + "Matrix": { + "Package": "Matrix", + "Version": "1.3-3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "df57c82e79600601287edfdcef92c2d6" + }, + "R6": { + "Package": "R6", + "Version": "2.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b203113193e70978a696b2809525649d" + }, + "RColorBrewer": { + "Package": "RColorBrewer", + "Version": "1.1-2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e031418365a7f7a766181ab5a41a5716" + }, + "Rcpp": { + "Package": "Rcpp", + "Version": "1.0.6", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "dbb5e436998a7eba5a9d682060533338" + }, + "RcppArmadillo": { + "Package": "RcppArmadillo", + "Version": "0.10.5.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b1c9574190c1c84417eedaa5901138bf" + }, + "RcppEigen": { + "Package": "RcppEigen", + "Version": "0.3.3.9.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ddfa72a87fdf4c80466a20818be91d00" + }, + "RcppParallel": { + "Package": "RcppParallel", + "Version": "5.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0325b5be38a02d192828027983c7b470" + }, + "StanHeaders": { + "Package": "StanHeaders", + "Version": "2.21.0-7", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0459d4dd7a8c239be18469a30c23dd4b" + }, + "V8": { + "Package": "V8", + "Version": "3.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b181d4eb2904da0897a742a559964173" + }, + "abind": { + "Package": "abind", + "Version": "1.4-5", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "4f57884290cc75ab22f4af9e9d4ca862" + }, + "askpass": { + "Package": "askpass", + "Version": "1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e8a22846fff485f0be3770c2da758713" + }, + "assertthat": { + "Package": "assertthat", + "Version": "0.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "50c838a310445e954bc13f26f26a6ecf" + }, + "backports": { + "Package": "backports", + "Version": "1.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "644043219fc24e190c2f620c1a380a69" + }, + "base64enc": { + "Package": "base64enc", + "Version": "0.1-3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "543776ae6848fde2f48ff3816d0628bc" + }, + "bayesplot": { + "Package": "bayesplot", + "Version": "1.8.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0a8b771b1b255fc55f271ae49e671afc" + }, + "blob": { + "Package": "blob", + "Version": "1.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "9addc7e2c5954eca5719928131fed98c" + }, + "boot": { + "Package": "boot", + "Version": "1.3-28", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0baa960e3b49c6176a4f42addcbacc59" + }, + "bridgesampling": { + "Package": "bridgesampling", + "Version": "1.1-2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a64d4be61a31696fec80dc24cef69a9c" + }, + "brms": { + "Package": "brms", + "Version": "2.15.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "fc74baf33468122b49024ba32e7c1a6a" + }, + "broom": { + "Package": "broom", + "Version": "0.7.8", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "c6dddf90a1e38665c211932696ddbfde" + }, + "bslib": { + "Package": "bslib", + "Version": "0.2.5.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2f069f3f42847231aef7baa49bed97b0" + }, + "cachem": { + "Package": "cachem", + "Version": "1.0.5", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "5346f76a33eb7417812c270b04a5581b" + }, + "callr": { + "Package": "callr", + "Version": "3.7.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "461aa75a11ce2400245190ef5d3995df" + }, + "cellranger": { + "Package": "cellranger", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "f61dbaec772ccd2e17705c1e872e9e7c" + }, + "checkmate": { + "Package": "checkmate", + "Version": "2.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a667800d5f0350371bedeb8b8b950289" + }, + "cli": { + "Package": "cli", + "Version": "2.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a94ba44cee3ea571e813721e64184172" + }, + "clipr": { + "Package": "clipr", + "Version": "0.7.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ebaa97ac99cc2daf04e77eecc7b781d7" + }, + "coda": { + "Package": "coda", + "Version": "0.19-4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "24b6d006b8b2343876cf230687546932" + }, + "codetools": { + "Package": "codetools", + "Version": "0.2-18", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "019388fc48e48b3da0d3a76ff94608a8" + }, + "colorspace": { + "Package": "colorspace", + "Version": "2.0-2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6baccb763ee83c0bd313460fdb8b8a84" + }, + "colourpicker": { + "Package": "colourpicker", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "fe0cb3d8854c168aef827aeab4b3b473" + }, + "commonmark": { + "Package": "commonmark", + "Version": "1.7", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0f22be39ec1d141fd03683c06f3a6e67" + }, + "cpp11": { + "Package": "cpp11", + "Version": "0.3.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e02edab2bc389c5e4b12949b13df44f2" + }, + "crayon": { + "Package": "crayon", + "Version": "1.4.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e75525c55c70e5f4f78c9960a4b402e9" + }, + "crosstalk": { + "Package": "crosstalk", + "Version": "1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2b06f9e415a62b6762e4b8098d2aecbc" + }, + "curl": { + "Package": "curl", + "Version": "4.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "022c42d49c28e95d69ca60446dbabf88" + }, + "data.table": { + "Package": "data.table", + "Version": "1.14.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d1b8b1a821ee564a3515fa6c6d5c52dc" + }, + "dbplyr": { + "Package": "dbplyr", + "Version": "2.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "1f37fa4ab2f5f7eded42f78b9a887182" + }, + "desc": { + "Package": "desc", + "Version": "1.3.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b6963166f7f10b970af1006c462ce6cd" + }, + "digest": { + "Package": "digest", + "Version": "0.6.27", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a0cbe758a531d054b537d16dff4d58a1" + }, + "dplyr": { + "Package": "dplyr", + "Version": "1.0.7", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "36f1ae62f026c8ba9f9b5c9a08c03297" + }, + "dtplyr": { + "Package": "dtplyr", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "1e14e4c5b2814de5225312394bc316da" + }, + "dygraphs": { + "Package": "dygraphs", + "Version": "1.1.1.6", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "716869fffc16e282c118f8894e082a7d" + }, + "ellipsis": { + "Package": "ellipsis", + "Version": "0.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "bb0eec2fe32e88d9e2836c2f73ea2077" + }, + "evaluate": { + "Package": "evaluate", + "Version": "0.14", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ec8ca05cffcc70569eaaad8469d2a3a7" + }, + "fansi": { + "Package": "fansi", + "Version": "0.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d447b40982c576a72b779f0a3b3da227" + }, + "farver": { + "Package": "farver", + "Version": "2.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "c98eb5133d9cb9e1622b8691487f11bb" + }, + "fastmap": { + "Package": "fastmap", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "77bd60a6157420d4ffa93b27cf6a58b8" + }, + "forcats": { + "Package": "forcats", + "Version": "0.5.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "81c3244cab67468aac4c60550832655d" + }, + "foreign": { + "Package": "foreign", + "Version": "0.8-81", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "74628ea7a3be5ee8a7b5bb0a8e84882e" + }, + "fs": { + "Package": "fs", + "Version": "1.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "44594a07a42e5f91fac9f93fda6d0109" + }, + "future": { + "Package": "future", + "Version": "1.21.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "f25fad6bee82b7ab01f055e2d813b96f" + }, + "gamm4": { + "Package": "gamm4", + "Version": "0.2-6", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "07fd04f71ff5fbf2a69c74e4ec128d65" + }, + "gargle": { + "Package": "gargle", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "bc3bc77ea458de0a6efe94e8e7e9c641" + }, + "generics": { + "Package": "generics", + "Version": "0.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "4d243a9c10b00589889fe32314ffd902" + }, + "ggplot2": { + "Package": "ggplot2", + "Version": "3.3.5", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d7566c471c7b17e095dd023b9ef155ad" + }, + "ggridges": { + "Package": "ggridges", + "Version": "0.5.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "9d028e8f37c84dba356ce3c367a1978e" + }, + "globals": { + "Package": "globals", + "Version": "0.14.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "eca8023ed5ca6372479ebb9b3207f5ae" + }, + "glue": { + "Package": "glue", + "Version": "1.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6efd734b14c6471cfe443345f3e35e29" + }, + "googledrive": { + "Package": "googledrive", + "Version": "1.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "79ba5d18133290a69b7c135dc3dfef1a" + }, + "googlesheets4": { + "Package": "googlesheets4", + "Version": "0.3.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "fd028aa45f6785d7e83c0ed896b45904" + }, + "gridExtra": { + "Package": "gridExtra", + "Version": "2.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "7d7f283939f563670a697165b2cf5560" + }, + "gtable": { + "Package": "gtable", + "Version": "0.3.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ac5c6baf7822ce8732b343f14c072c4d" + }, + "gtools": { + "Package": "gtools", + "Version": "3.9.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2ace6c4a06297d0b364e0444384a2b82" + }, + "haven": { + "Package": "haven", + "Version": "2.4.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "9ace652bffef34af558eb5d70903ac2f" + }, + "highr": { + "Package": "highr", + "Version": "0.9", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "8eb36c8125038e648e5d111c0d7b2ed4" + }, + "hms": { + "Package": "hms", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e4bf161ccb74a2c1c0e8ac63bbe332b4" + }, + "htmltools": { + "Package": "htmltools", + "Version": "0.5.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "af2c2531e55df5cf230c4b5444fc973c" + }, + "htmlwidgets": { + "Package": "htmlwidgets", + "Version": "1.5.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6fdaa86d0700f8b3e92ee3c445a5a10d" + }, + "httpuv": { + "Package": "httpuv", + "Version": "1.6.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "54344a78aae37bc6ef39b1240969df8e" + }, + "httr": { + "Package": "httr", + "Version": "1.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a525aba14184fec243f9eaec62fbed43" + }, + "ids": { + "Package": "ids", + "Version": "1.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "99df65cfef20e525ed38c3d2577f7190" + }, + "igraph": { + "Package": "igraph", + "Version": "1.2.6", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "7b1f856410253d56ea67ad808f7cdff6" + }, + "inline": { + "Package": "inline", + "Version": "0.3.19", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "1deaf1de3eac7e1d3377954b3a283652" + }, + "isoband": { + "Package": "isoband", + "Version": "0.2.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b2008df40fb297e3fef135c7e8eeec1a" + }, + "jquerylib": { + "Package": "jquerylib", + "Version": "0.1.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "5aab57a3bd297eee1c1d862735972182" + }, + "jsonlite": { + "Package": "jsonlite", + "Version": "1.7.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "98138e0994d41508c7a6b84a0600cfcb" + }, + "knitr": { + "Package": "knitr", + "Version": "1.33", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0bc1b5da1b0eb07cd4b727e95e9ff0b8" + }, + "labeling": { + "Package": "labeling", + "Version": "0.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "3d5108641f47470611a32d0bdf357a72" + }, + "later": { + "Package": "later", + "Version": "1.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b61890ae77fea19fc8acadd25db70aa4" + }, + "lattice": { + "Package": "lattice", + "Version": "0.20-44", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "f36bf1a849d9106dc2af72e501f9de41" + }, + "lazyeval": { + "Package": "lazyeval", + "Version": "0.2.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d908914ae53b04d4c0c0fd72ecc35370" + }, + "lifecycle": { + "Package": "lifecycle", + "Version": "1.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "3471fb65971f1a7b2d4ae7848cf2db8d" + }, + "listenv": { + "Package": "listenv", + "Version": "0.8.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0bde42ee282efb18c7c4e63822f5b4f7" + }, + "lme4": { + "Package": "lme4", + "Version": "1.1-27.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "c995b0405ce0894d6fe52b3e08ea9085" + }, + "loo": { + "Package": "loo", + "Version": "2.4.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "cb08bd08c899e7d68661714d699c42e4" + }, + "lubridate": { + "Package": "lubridate", + "Version": "1.7.10", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "1ebfdc8a3cfe8fe19184f5481972b092" + }, + "magrittr": { + "Package": "magrittr", + "Version": "2.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "41287f1ac7d28a92f0a286ed507928d3" + }, + "markdown": { + "Package": "markdown", + "Version": "1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "61e4a10781dd00d7d81dd06ca9b94e95" + }, + "matrixStats": { + "Package": "matrixStats", + "Version": "0.59.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "45a2fee40a2f1bcad392cbe8a5f10d4f" + }, + "mgcv": { + "Package": "mgcv", + "Version": "1.8-35", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "89fd8b2ad4a6cb4979b78cf2a77ab503" + }, + "mime": { + "Package": "mime", + "Version": "0.11", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "8974a907200fc9948d636fe7d85ca9fb" + }, + "miniUI": { + "Package": "miniUI", + "Version": "0.1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "fec5f52652d60615fdb3957b3d74324a" + }, + "minqa": { + "Package": "minqa", + "Version": "1.2.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "eaee7d2a6f3ed4491df868611cb064cc" + }, + "modelr": { + "Package": "modelr", + "Version": "0.1.8", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "9fd59716311ee82cba83dc2826fc5577" + }, + "munsell": { + "Package": "munsell", + "Version": "0.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6dfe8bf774944bd5595785e3229d8771" + }, + "mvtnorm": { + "Package": "mvtnorm", + "Version": "1.1-2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6f0133c3842aef0394dbd844a21d3f5f" + }, + "nleqslv": { + "Package": "nleqslv", + "Version": "3.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "bd00b0d9f6e1aa1462a89eb05d992b38" + }, + "nlme": { + "Package": "nlme", + "Version": "3.1-152", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "35de1ce639f20b5e10f7f46260730c65" + }, + "nloptr": { + "Package": "nloptr", + "Version": "1.2.2.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2737faeee353704efec5afa1e943dd64" + }, + "numDeriv": { + "Package": "numDeriv", + "Version": "2016.8-1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "df58958f293b166e4ab885ebcad90e02" + }, + "openssl": { + "Package": "openssl", + "Version": "1.4.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "f4dbc5a47fd93d3415249884d31d6791" + }, + "optimx": { + "Package": "optimx", + "Version": "2021-6.12", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "dee57179a02ea2fcdf0aa32a9f2114b7" + }, + "packrat": { + "Package": "packrat", + "Version": "0.6.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0d6cc4c357e7602bb3eee299f4cfc2a5" + }, + "parallelly": { + "Package": "parallelly", + "Version": "1.26.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "08925887608f7086dd7522933ff43a93" + }, + "pillar": { + "Package": "pillar", + "Version": "1.6.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "8672ae02bd20f6479bce2d06c7ff1401" + }, + "pkgbuild": { + "Package": "pkgbuild", + "Version": "1.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "725fcc30222d4d11ec68efb8ff11a9af" + }, + "pkgconfig": { + "Package": "pkgconfig", + "Version": "2.0.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "01f28d4278f15c76cddbea05899c5d6f" + }, + "plyr": { + "Package": "plyr", + "Version": "1.8.6", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ec0e5ab4e5f851f6ef32cd1d1984957f" + }, + "prettyunits": { + "Package": "prettyunits", + "Version": "1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "95ef9167b75dde9d2ccc3c7528393e7e" + }, + "processx": { + "Package": "processx", + "Version": "3.5.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0cbca2bc4d16525d009c4dbba156b37c" + }, + "progress": { + "Package": "progress", + "Version": "1.2.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "14dc9f7a3c91ebb14ec5bb9208a07061" + }, + "projpred": { + "Package": "projpred", + "Version": "2.0.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "69ee4853dcc10a5966905c26c6e0fbab" + }, + "promises": { + "Package": "promises", + "Version": "1.2.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "4ab2c43adb4d4699cf3690acd378d75d" + }, + "ps": { + "Package": "ps", + "Version": "1.6.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "32620e2001c1dce1af49c49dccbb9420" + }, + "purrr": { + "Package": "purrr", + "Version": "0.3.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "97def703420c8ab10d8f0e6c72101e02" + }, + "rappdirs": { + "Package": "rappdirs", + "Version": "0.3.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "5e3c5dc0b071b21fa128676560dbe94d" + }, + "readr": { + "Package": "readr", + "Version": "1.4.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2639976851f71f330264a9c9c3d43a61" + }, + "readxl": { + "Package": "readxl", + "Version": "1.3.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "63537c483c2dbec8d9e3183b3735254a" + }, + "rematch": { + "Package": "rematch", + "Version": "1.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "c66b930d20bb6d858cd18e1cebcfae5c" + }, + "rematch2": { + "Package": "rematch2", + "Version": "2.1.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "76c9e04c712a05848ae7a23d2f170a40" + }, + "renv": { + "Package": "renv", + "Version": "0.13.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "079cb1f03ff972b30401ed05623cbe92" + }, + "reprex": { + "Package": "reprex", + "Version": "2.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "8482bbeef0c194ac236aef7c51ee375f" + }, + "reshape2": { + "Package": "reshape2", + "Version": "1.4.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "bb5996d0bd962d214a11140d77589917" + }, + "rlang": { + "Package": "rlang", + "Version": "0.4.11", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "515f341d3affe0de9e4a7f762efb0456" + }, + "rmarkdown": { + "Package": "rmarkdown", + "Version": "2.9", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "912c09266d5470516df4df7a303cde92" + }, + "rngtools": { + "Package": "rngtools", + "Version": "1.5", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ee3ca1cacd13e35dbb71ac1999d98cbc" + }, + "rprojroot": { + "Package": "rprojroot", + "Version": "2.0.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "249d8cd1e74a8f6a26194a91b47f21d1" + }, + "rsconnect": { + "Package": "rsconnect", + "Version": "0.8.18", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "50b4e52011e335f7d3053cfde69893ba" + }, + "rstan": { + "Package": "rstan", + "Version": "2.21.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "52772d81aa532a6331fd535701882c12" + }, + "rstantools": { + "Package": "rstantools", + "Version": "2.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "c91e0f20c967e246cb6f2efe8c60e15b" + }, + "rstudioapi": { + "Package": "rstudioapi", + "Version": "0.13", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "06c85365a03fdaf699966cc1d3cf53ea" + }, + "rvest": { + "Package": "rvest", + "Version": "1.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "74b905b0076e1de6e27f540c95ba68d5" + }, + "sass": { + "Package": "sass", + "Version": "0.4.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "50cf822feb64bb3977bda0b7091be623" + }, + "scales": { + "Package": "scales", + "Version": "1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6f76f71042411426ec8df6c54f34e6dd" + }, + "selectr": { + "Package": "selectr", + "Version": "0.4-2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "3838071b66e0c566d55cc26bd6e27bf4" + }, + "shiny": { + "Package": "shiny", + "Version": "1.6.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "6e3b6ae7fe02b5859e4bb277f218b8ae" + }, + "shinyjs": { + "Package": "shinyjs", + "Version": "2.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "9ddfc91d4280eaa34c2103951538976f" + }, + "shinystan": { + "Package": "shinystan", + "Version": "2.5.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "85e1a9e77cdd1b8740e92b37f8fdce7b" + }, + "shinythemes": { + "Package": "shinythemes", + "Version": "1.2.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "30f0ebc41feba25691073626ff5e2cf4" + }, + "sourcetools": { + "Package": "sourcetools", + "Version": "0.1.7", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "947e4e02a79effa5d512473e10f41797" + }, + "stringi": { + "Package": "stringi", + "Version": "1.6.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "9df5e6f9a7fa11b84adf0429961de66a" + }, + "stringr": { + "Package": "stringr", + "Version": "1.4.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "0759e6b6c0957edb1311028a49a35e76" + }, + "sys": { + "Package": "sys", + "Version": "3.4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b227d13e29222b4574486cfcbde077fa" + }, + "threejs": { + "Package": "threejs", + "Version": "0.3.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2ad32c3a8745e827977f394bc387e3b0" + }, + "tibble": { + "Package": "tibble", + "Version": "3.1.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "349b40a9f144516d537c875e786ec8b8" + }, + "tidyr": { + "Package": "tidyr", + "Version": "1.1.3", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "450d7dfaedde58e28586b854eeece4fa" + }, + "tidyselect": { + "Package": "tidyselect", + "Version": "1.1.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "7243004a708d06d4716717fa1ff5b2fe" + }, + "tidyverse": { + "Package": "tidyverse", + "Version": "1.3.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "fc4c72b6ae9bb283416bd59a3303bbab" + }, + "tinytex": { + "Package": "tinytex", + "Version": "0.32", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "db9a6f2cf147751322d22c9f6647c7bd" + }, + "utf8": { + "Package": "utf8", + "Version": "1.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "c3ad47dc6da0751f18ed53c4613e3ac7" + }, + "uuid": { + "Package": "uuid", + "Version": "0.1-4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "e4169eb989a5d03ccb6b628cad1b1b50" + }, + "vctrs": { + "Package": "vctrs", + "Version": "0.3.8", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ecf749a1b39ea72bd9b51b76292261f1" + }, + "viridisLite": { + "Package": "viridisLite", + "Version": "0.4.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "55e157e2aa88161bdb0754218470d204" + }, + "withr": { + "Package": "withr", + "Version": "2.4.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ad03909b44677f930fa156d47d7a3aeb" + }, + "xfun": { + "Package": "xfun", + "Version": "0.24", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "88cdb9779a657ad80ad942245fffba31" + }, + "xml2": { + "Package": "xml2", + "Version": "1.3.2", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "d4d71a75dd3ea9eb5fa28cc21f9585e2" + }, + "xtable": { + "Package": "xtable", + "Version": "1.8-4", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "b8acdf8af494d9ec19ccb2481a9b11c2" + }, + "xts": { + "Package": "xts", + "Version": "0.12.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "ca2fd4ad8ef78cca3aa2b30f992798a8" + }, + "yaml": { + "Package": "yaml", + "Version": "2.2.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "2826c5d9efb0a88f657c7a679c7106db" + }, + "zoo": { + "Package": "zoo", + "Version": "1.8-9", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "035d1c7c12593038c26fb1c2fd40c4d2" + } + } +} diff --git a/renv/.gitignore b/renv/.gitignore new file mode 100644 index 0000000..2129631 --- /dev/null +++ b/renv/.gitignore @@ -0,0 +1,5 @@ +library/ +local/ +lock/ +python/ +staging/ diff --git a/renv/activate.R b/renv/activate.R new file mode 100644 index 0000000..b852628 --- /dev/null +++ b/renv/activate.R @@ -0,0 +1,654 @@ + +local({ + + # the requested version of renv + version <- "0.13.2" + + # the project directory + project <- getwd() + + # avoid recursion + if (!is.na(Sys.getenv("RENV_R_INITIALIZING", unset = NA))) + return(invisible(TRUE)) + + # signal that we're loading renv during R startup + Sys.setenv("RENV_R_INITIALIZING" = "true") + on.exit(Sys.unsetenv("RENV_R_INITIALIZING"), add = TRUE) + + # signal that we've consented to use renv + options(renv.consent = TRUE) + + # load the 'utils' package eagerly -- this ensures that renv shims, which + # mask 'utils' packages, will come first on the search path + library(utils, lib.loc = .Library) + + # check to see if renv has already been loaded + if ("renv" %in% loadedNamespaces()) { + + # if renv has already been loaded, and it's the requested version of renv, + # nothing to do + spec <- .getNamespaceInfo(.getNamespace("renv"), "spec") + if (identical(spec[["version"]], version)) + return(invisible(TRUE)) + + # otherwise, unload and attempt to load the correct version of renv + unloadNamespace("renv") + + } + + # load bootstrap tools + bootstrap <- function(version, library) { + + # attempt to download renv + tarball <- tryCatch(renv_bootstrap_download(version), error = identity) + if (inherits(tarball, "error")) + stop("failed to download renv ", version) + + # now attempt to install + status <- tryCatch(renv_bootstrap_install(version, tarball, library), error = identity) + if (inherits(status, "error")) + stop("failed to install renv ", version) + + } + + renv_bootstrap_tests_running <- function() { + getOption("renv.tests.running", default = FALSE) + } + + renv_bootstrap_repos <- function() { + + # check for repos override + repos <- Sys.getenv("RENV_CONFIG_REPOS_OVERRIDE", unset = NA) + if (!is.na(repos)) + return(repos) + + # if we're testing, re-use the test repositories + if (renv_bootstrap_tests_running()) + return(getOption("renv.tests.repos")) + + # retrieve current repos + repos <- getOption("repos") + + # ensure @CRAN@ entries are resolved + repos[repos == "@CRAN@"] <- getOption( + "renv.repos.cran", + "https://cloud.r-project.org" + ) + + # add in renv.bootstrap.repos if set + default <- c(FALLBACK = "https://cloud.r-project.org") + extra <- getOption("renv.bootstrap.repos", default = default) + repos <- c(repos, extra) + + # remove duplicates that might've snuck in + dupes <- duplicated(repos) | duplicated(names(repos)) + repos[!dupes] + + } + + renv_bootstrap_download <- function(version) { + + # if the renv version number has 4 components, assume it must + # be retrieved via github + nv <- numeric_version(version) + components <- unclass(nv)[[1]] + + methods <- if (length(components) == 4L) { + list( + renv_bootstrap_download_github + ) + } else { + list( + renv_bootstrap_download_cran_latest, + renv_bootstrap_download_cran_archive + ) + } + + for (method in methods) { + path <- tryCatch(method(version), error = identity) + if (is.character(path) && file.exists(path)) + return(path) + } + + stop("failed to download renv ", version) + + } + + renv_bootstrap_download_impl <- function(url, destfile) { + + mode <- "wb" + + # https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17715 + fixup <- + Sys.info()[["sysname"]] == "Windows" && + substring(url, 1L, 5L) == "file:" + + if (fixup) + mode <- "w+b" + + utils::download.file( + url = url, + destfile = destfile, + mode = mode, + quiet = TRUE + ) + + } + + renv_bootstrap_download_cran_latest <- function(version) { + + spec <- renv_bootstrap_download_cran_latest_find(version) + + message("* Downloading renv ", version, " ... ", appendLF = FALSE) + + type <- spec$type + repos <- spec$repos + + info <- tryCatch( + utils::download.packages( + pkgs = "renv", + destdir = tempdir(), + repos = repos, + type = type, + quiet = TRUE + ), + condition = identity + ) + + if (inherits(info, "condition")) { + message("FAILED") + return(FALSE) + } + + # report success and return + message("OK (downloaded ", type, ")") + info[1, 2] + + } + + renv_bootstrap_download_cran_latest_find <- function(version) { + + # check whether binaries are supported on this system + binary <- + getOption("renv.bootstrap.binary", default = TRUE) && + !identical(.Platform$pkgType, "source") && + !identical(getOption("pkgType"), "source") && + Sys.info()[["sysname"]] %in% c("Darwin", "Windows") + + types <- c(if (binary) "binary", "source") + + # iterate over types + repositories + for (type in types) { + for (repos in renv_bootstrap_repos()) { + + # retrieve package database + db <- tryCatch( + as.data.frame( + utils::available.packages(type = type, repos = repos), + stringsAsFactors = FALSE + ), + error = identity + ) + + if (inherits(db, "error")) + next + + # check for compatible entry + entry <- db[db$Package %in% "renv" & db$Version %in% version, ] + if (nrow(entry) == 0) + next + + # found it; return spec to caller + spec <- list(entry = entry, type = type, repos = repos) + return(spec) + + } + } + + # if we got here, we failed to find renv + fmt <- "renv %s is not available from your declared package repositories" + stop(sprintf(fmt, version)) + + } + + renv_bootstrap_download_cran_archive <- function(version) { + + name <- sprintf("renv_%s.tar.gz", version) + repos <- renv_bootstrap_repos() + urls <- file.path(repos, "src/contrib/Archive/renv", name) + destfile <- file.path(tempdir(), name) + + message("* Downloading renv ", version, " ... ", appendLF = FALSE) + + for (url in urls) { + + status <- tryCatch( + renv_bootstrap_download_impl(url, destfile), + condition = identity + ) + + if (identical(status, 0L)) { + message("OK") + return(destfile) + } + + } + + message("FAILED") + return(FALSE) + + } + + renv_bootstrap_download_github <- function(version) { + + enabled <- Sys.getenv("RENV_BOOTSTRAP_FROM_GITHUB", unset = "TRUE") + if (!identical(enabled, "TRUE")) + return(FALSE) + + # prepare download options + pat <- Sys.getenv("GITHUB_PAT") + if (nzchar(Sys.which("curl")) && nzchar(pat)) { + fmt <- "--location --fail --header \"Authorization: token %s\"" + extra <- sprintf(fmt, pat) + saved <- options("download.file.method", "download.file.extra") + options(download.file.method = "curl", download.file.extra = extra) + on.exit(do.call(base::options, saved), add = TRUE) + } else if (nzchar(Sys.which("wget")) && nzchar(pat)) { + fmt <- "--header=\"Authorization: token %s\"" + extra <- sprintf(fmt, pat) + saved <- options("download.file.method", "download.file.extra") + options(download.file.method = "wget", download.file.extra = extra) + on.exit(do.call(base::options, saved), add = TRUE) + } + + message("* Downloading renv ", version, " from GitHub ... ", appendLF = FALSE) + + url <- file.path("https://api.github.com/repos/rstudio/renv/tarball", version) + name <- sprintf("renv_%s.tar.gz", version) + destfile <- file.path(tempdir(), name) + + status <- tryCatch( + renv_bootstrap_download_impl(url, destfile), + condition = identity + ) + + if (!identical(status, 0L)) { + message("FAILED") + return(FALSE) + } + + message("OK") + return(destfile) + + } + + renv_bootstrap_install <- function(version, tarball, library) { + + # attempt to install it into project library + message("* Installing renv ", version, " ... ", appendLF = FALSE) + dir.create(library, showWarnings = FALSE, recursive = TRUE) + + # invoke using system2 so we can capture and report output + bin <- R.home("bin") + exe <- if (Sys.info()[["sysname"]] == "Windows") "R.exe" else "R" + r <- file.path(bin, exe) + args <- c("--vanilla", "CMD", "INSTALL", "-l", shQuote(library), shQuote(tarball)) + output <- system2(r, args, stdout = TRUE, stderr = TRUE) + message("Done!") + + # check for successful install + status <- attr(output, "status") + if (is.numeric(status) && !identical(status, 0L)) { + header <- "Error installing renv:" + lines <- paste(rep.int("=", nchar(header)), collapse = "") + text <- c(header, lines, output) + writeLines(text, con = stderr()) + } + + status + + } + + renv_bootstrap_platform_prefix <- function() { + + # construct version prefix + version <- paste(R.version$major, R.version$minor, sep = ".") + prefix <- paste("R", numeric_version(version)[1, 1:2], sep = "-") + + # include SVN revision for development versions of R + # (to avoid sharing platform-specific artefacts with released versions of R) + devel <- + identical(R.version[["status"]], "Under development (unstable)") || + identical(R.version[["nickname"]], "Unsuffered Consequences") + + if (devel) + prefix <- paste(prefix, R.version[["svn rev"]], sep = "-r") + + # build list of path components + components <- c(prefix, R.version$platform) + + # include prefix if provided by user + prefix <- renv_bootstrap_platform_prefix_impl() + if (!is.na(prefix) && nzchar(prefix)) + components <- c(prefix, components) + + # build prefix + paste(components, collapse = "/") + + } + + renv_bootstrap_platform_prefix_impl <- function() { + + # if an explicit prefix has been supplied, use it + prefix <- Sys.getenv("RENV_PATHS_PREFIX", unset = NA) + if (!is.na(prefix)) + return(prefix) + + # if the user has requested an automatic prefix, generate it + auto <- Sys.getenv("RENV_PATHS_PREFIX_AUTO", unset = NA) + if (auto %in% c("TRUE", "True", "true", "1")) + return(renv_bootstrap_platform_prefix_auto()) + + # empty string on failure + "" + + } + + renv_bootstrap_platform_prefix_auto <- function() { + + prefix <- tryCatch(renv_bootstrap_platform_os(), error = identity) + if (inherits(prefix, "error") || prefix %in% "unknown") { + + msg <- paste( + "failed to infer current operating system", + "please file a bug report at https://github.com/rstudio/renv/issues", + sep = "; " + ) + + warning(msg) + + } + + prefix + + } + + renv_bootstrap_platform_os <- function() { + + sysinfo <- Sys.info() + sysname <- sysinfo[["sysname"]] + + # handle Windows + macOS up front + if (sysname == "Windows") + return("windows") + else if (sysname == "Darwin") + return("macos") + + # check for os-release files + for (file in c("/etc/os-release", "/usr/lib/os-release")) + if (file.exists(file)) + return(renv_bootstrap_platform_os_via_os_release(file, sysinfo)) + + # check for redhat-release files + if (file.exists("/etc/redhat-release")) + return(renv_bootstrap_platform_os_via_redhat_release()) + + "unknown" + + } + + renv_bootstrap_platform_os_via_os_release <- function(file, sysinfo) { + + # read /etc/os-release + release <- utils::read.table( + file = file, + sep = "=", + quote = c("\"", "'"), + col.names = c("Key", "Value"), + comment.char = "#", + stringsAsFactors = FALSE + ) + + vars <- as.list(release$Value) + names(vars) <- release$Key + + # get os name + os <- tolower(sysinfo[["sysname"]]) + + # read id + id <- "unknown" + for (field in c("ID", "ID_LIKE")) { + if (field %in% names(vars) && nzchar(vars[[field]])) { + id <- vars[[field]] + break + } + } + + # read version + version <- "unknown" + for (field in c("UBUNTU_CODENAME", "VERSION_CODENAME", "VERSION_ID", "BUILD_ID")) { + if (field %in% names(vars) && nzchar(vars[[field]])) { + version <- vars[[field]] + break + } + } + + # join together + paste(c(os, id, version), collapse = "-") + + } + + renv_bootstrap_platform_os_via_redhat_release <- function() { + + # read /etc/redhat-release + contents <- readLines("/etc/redhat-release", warn = FALSE) + + # infer id + id <- if (grepl("centos", contents, ignore.case = TRUE)) + "centos" + else if (grepl("redhat", contents, ignore.case = TRUE)) + "redhat" + else + "unknown" + + # try to find a version component (very hacky) + version <- "unknown" + + parts <- strsplit(contents, "[[:space:]]")[[1L]] + for (part in parts) { + + nv <- tryCatch(numeric_version(part), error = identity) + if (inherits(nv, "error")) + next + + version <- nv[1, 1] + break + + } + + paste(c("linux", id, version), collapse = "-") + + } + + renv_bootstrap_library_root_name <- function(project) { + + # use project name as-is if requested + asis <- Sys.getenv("RENV_PATHS_LIBRARY_ROOT_ASIS", unset = "FALSE") + if (asis) + return(basename(project)) + + # otherwise, disambiguate based on project's path + id <- substring(renv_bootstrap_hash_text(project), 1L, 8L) + paste(basename(project), id, sep = "-") + + } + + renv_bootstrap_library_root <- function(project) { + + path <- Sys.getenv("RENV_PATHS_LIBRARY", unset = NA) + if (!is.na(path)) + return(path) + + path <- Sys.getenv("RENV_PATHS_LIBRARY_ROOT", unset = NA) + if (!is.na(path)) { + name <- renv_bootstrap_library_root_name(project) + return(file.path(path, name)) + } + + prefix <- renv_bootstrap_profile_prefix() + paste(c(project, prefix, "renv/library"), collapse = "/") + + } + + renv_bootstrap_validate_version <- function(version) { + + loadedversion <- utils::packageDescription("renv", fields = "Version") + if (version == loadedversion) + return(TRUE) + + # assume four-component versions are from GitHub; three-component + # versions are from CRAN + components <- strsplit(loadedversion, "[.-]")[[1]] + remote <- if (length(components) == 4L) + paste("rstudio/renv", loadedversion, sep = "@") + else + paste("renv", loadedversion, sep = "@") + + fmt <- paste( + "renv %1$s was loaded from project library, but this project is configured to use renv %2$s.", + "Use `renv::record(\"%3$s\")` to record renv %1$s in the lockfile.", + "Use `renv::restore(packages = \"renv\")` to install renv %2$s into the project library.", + sep = "\n" + ) + + msg <- sprintf(fmt, loadedversion, version, remote) + warning(msg, call. = FALSE) + + FALSE + + } + + renv_bootstrap_hash_text <- function(text) { + + hashfile <- tempfile("renv-hash-") + on.exit(unlink(hashfile), add = TRUE) + + writeLines(text, con = hashfile) + tools::md5sum(hashfile) + + } + + renv_bootstrap_load <- function(project, libpath, version) { + + # try to load renv from the project library + if (!requireNamespace("renv", lib.loc = libpath, quietly = TRUE)) + return(FALSE) + + # warn if the version of renv loaded does not match + renv_bootstrap_validate_version(version) + + # load the project + renv::load(project) + + TRUE + + } + + renv_bootstrap_profile_load <- function(project) { + + # if RENV_PROFILE is already set, just use that + profile <- Sys.getenv("RENV_PROFILE", unset = NA) + if (!is.na(profile) && nzchar(profile)) + return(profile) + + # check for a profile file (nothing to do if it doesn't exist) + path <- file.path(project, "renv/local/profile") + if (!file.exists(path)) + return(NULL) + + # read the profile, and set it if it exists + contents <- readLines(path, warn = FALSE) + if (length(contents) == 0L) + return(NULL) + + # set RENV_PROFILE + profile <- contents[[1L]] + if (nzchar(profile)) + Sys.setenv(RENV_PROFILE = profile) + + profile + + } + + renv_bootstrap_profile_prefix <- function() { + profile <- renv_bootstrap_profile_get() + if (!is.null(profile)) + return(file.path("renv/profiles", profile)) + } + + renv_bootstrap_profile_get <- function() { + profile <- Sys.getenv("RENV_PROFILE", unset = "") + renv_bootstrap_profile_normalize(profile) + } + + renv_bootstrap_profile_set <- function(profile) { + profile <- renv_bootstrap_profile_normalize(profile) + if (is.null(profile)) + Sys.unsetenv("RENV_PROFILE") + else + Sys.setenv(RENV_PROFILE = profile) + } + + renv_bootstrap_profile_normalize <- function(profile) { + + if (is.null(profile) || profile %in% c("", "default")) + return(NULL) + + profile + + } + + # load the renv profile, if any + renv_bootstrap_profile_load(project) + + # construct path to library root + root <- renv_bootstrap_library_root(project) + + # construct library prefix for platform + prefix <- renv_bootstrap_platform_prefix() + + # construct full libpath + libpath <- file.path(root, prefix) + + # attempt to load + if (renv_bootstrap_load(project, libpath, version)) + return(TRUE) + + # load failed; inform user we're about to bootstrap + prefix <- paste("# Bootstrapping renv", version) + postfix <- paste(rep.int("-", 77L - nchar(prefix)), collapse = "") + header <- paste(prefix, postfix) + message(header) + + # perform bootstrap + bootstrap(version, libpath) + + # exit early if we're just testing bootstrap + if (!is.na(Sys.getenv("RENV_BOOTSTRAP_INSTALL_ONLY", unset = NA))) + return(TRUE) + + # try again to load + if (requireNamespace("renv", lib.loc = libpath, quietly = TRUE)) { + message("* Successfully installed and loaded renv ", version, ".") + return(renv::load()) + } + + # failed to download or load renv; warn the user + msg <- c( + "Failed to find an renv installation: the project will not be loaded.", + "Use `renv::activate()` to re-initialize the project." + ) + + warning(paste(msg, collapse = "\n"), call. = FALSE) + +}) diff --git a/renv/settings.dcf b/renv/settings.dcf new file mode 100644 index 0000000..fc4e479 --- /dev/null +++ b/renv/settings.dcf @@ -0,0 +1,8 @@ +external.libraries: +ignored.packages: +package.dependency.fields: Imports, Depends, LinkingTo +r.version: +snapshot.type: implicit +use.cache: TRUE +vcs.ignore.library: TRUE +vcs.ignore.local: TRUE