Binary files /tmp/tmpah1mjzoy/8NqVPGu0hI/r-cran-mcmcpack-1.5-0/build/partial.rdb and /tmp/tmpah1mjzoy/VVx9PZA7pZ/r-cran-mcmcpack-1.6-0/build/partial.rdb differ Binary files /tmp/tmpah1mjzoy/8NqVPGu0hI/r-cran-mcmcpack-1.5-0/data/Euro2016.rda and /tmp/tmpah1mjzoy/VVx9PZA7pZ/r-cran-mcmcpack-1.6-0/data/Euro2016.rda differ diff -Nru r-cran-mcmcpack-1.5-0/debian/changelog r-cran-mcmcpack-1.6-0/debian/changelog --- r-cran-mcmcpack-1.5-0/debian/changelog 2021-01-24 20:10:32.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/debian/changelog 2021-10-14 20:24:18.000000000 +0000 @@ -1,3 +1,11 @@ +r-cran-mcmcpack (1.6-0-1) unstable; urgency=medium + + * Team upload. + * New upstream version + * Standards-Version: 4.6.0 (routine-update) + + -- Andreas Tille Thu, 14 Oct 2021 22:24:18 +0200 + r-cran-mcmcpack (1.5-0-1) unstable; urgency=medium * Team upload. diff -Nru r-cran-mcmcpack-1.5-0/debian/control r-cran-mcmcpack-1.6-0/debian/control --- r-cran-mcmcpack-1.5-0/debian/control 2021-01-24 20:10:32.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/debian/control 2021-10-14 20:24:18.000000000 +0000 @@ -12,7 +12,7 @@ r-cran-lattice, r-cran-mcmc, r-cran-quantreg -Standards-Version: 4.5.1 +Standards-Version: 4.6.0 Vcs-Browser: https://salsa.debian.org/r-pkg-team/r-cran-mcmcpack Vcs-Git: https://salsa.debian.org/r-pkg-team/r-cran-mcmcpack.git Homepage: https://cran.r-project.org/package=MCMCpack diff -Nru r-cran-mcmcpack-1.5-0/DESCRIPTION r-cran-mcmcpack-1.6-0/DESCRIPTION --- r-cran-mcmcpack-1.5-0/DESCRIPTION 2021-01-20 11:50:12.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/DESCRIPTION 2021-10-06 05:40:02.000000000 +0000 @@ -1,8 +1,8 @@ Package: MCMCpack -Version: 1.5-0 -Date: 2021-01-19 +Version: 1.6-0 +Date: 2021-10-05 Title: Markov Chain Monte Carlo (MCMC) Package -Author: Andrew D. Martin [aut], Kevin M. Quinn [aut], Jong Hee Park [aut,cre], Ghislain Vieilledent [ctb], Michael Malecki[ctb], Matthew Blackwell [ctb], Keith Poole [ctb], Craig Reed [ctb], Ben Goodrich [ctb], Ross Ihaka [cph], The R Development Core Team [cph], The R Foundation [cph], Pierre L'Ecuyer [cph], Makoto Matsumoto [cph], Takuji Nishimura [cph] +Author: Andrew D. Martin [aut], Kevin M. Quinn [aut], Jong Hee Park [aut,cre], Ghislain Vieilledent [ctb], Michael Malecki[ctb], Matthew Blackwell [ctb], Keith Poole [ctb], Craig Reed [ctb], Ben Goodrich [ctb], Qiushi Yu [ctb], Ross Ihaka [cph], The R Development Core Team [cph], The R Foundation [cph], Pierre L'Ecuyer [cph], Makoto Matsumoto [cph], Takuji Nishimura [cph] Maintainer: Jong Hee Park Depends: R (>= 3.6), coda (>= 0.11-3), MASS, stats Imports: graphics, grDevices, lattice, methods, utils, mcmc, quantreg @@ -19,9 +19,9 @@ License: GPL-3 SystemRequirements: gcc (>= 4.0) URL: https://CRAN.R-project.org/package=MCMCpack -Packaged: 2021-01-20 00:41:13 UTC; park +Packaged: 2021-10-06 02:04:17 UTC; park NeedsCompilation: yes Repository: CRAN RoxygenNote: 7.1.1 Encoding: UTF-8 -Date/Publication: 2021-01-20 11:50:12 UTC +Date/Publication: 2021-10-06 05:40:02 UTC diff -Nru r-cran-mcmcpack-1.5-0/inst/CITATION r-cran-mcmcpack-1.6-0/inst/CITATION --- r-cran-mcmcpack-1.5-0/inst/CITATION 2021-01-19 02:32:23.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/inst/CITATION 2021-10-06 01:35:27.000000000 +0000 @@ -10,12 +10,12 @@ volume = {42}, number = {9}, pages = {1--21}, - url = {"https://www.jstatsoft.org/v42/i09/"}, + Doi = {"10.18637/jss.v042.i09"}, textVersion = paste("Andrew D. Martin, Kevin M. Quinn, Jong Hee Park (2011).", "MCMCpack: Markov Chain Monte Carlo in R.", "Journal of Statistical Software. 42(9): 1-21.", - "URL https://www.jstatsoft.org/v42/i09/.") + "DOI 10.18637/jss.v042.i09.") ) diff -Nru r-cran-mcmcpack-1.5-0/inst/HISTORY.txt r-cran-mcmcpack-1.6-0/inst/HISTORY.txt --- r-cran-mcmcpack-1.5-0/inst/HISTORY.txt 2020-02-13 05:58:41.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/inst/HISTORY.txt 2021-07-25 01:56:03.000000000 +0000 @@ -1,6 +1,20 @@ // // Changes and Bug Fixes // + +1.5-0 to 1.6-0 + * MCMCpaircompare.R is added. + +1.4-9 to 1.5-0 + * Removed print(, digits=2) in MCMCmixfactanal + +1.4-8 to 1.4-9 + * error in marginal likelihood computation in MCMCregressChange + Thanks to Sid Chib for letting us know it + +1.4-6 to 1.4-7 + * print(, digits=2) in MCMCmixfactanal + 1.4-5 to 1.4-6 * Fixed ... in testSubjectBreak.Rd diff -Nru r-cran-mcmcpack-1.5-0/man/Euro2016.Rd r-cran-mcmcpack-1.6-0/man/Euro2016.Rd --- r-cran-mcmcpack-1.5-0/man/Euro2016.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/Euro2016.Rd 2021-07-25 01:53:29.000000000 +0000 @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{Euro2016} +\alias{Euro2016} +\title{Euro 2016 data} +\format{ +This dataframe contains all of the head-to-head results from +Euro 2016. This includes results from both the group stage and the +knock-out rounds. +\describe{ + \item{dummy.rater}{An artificial "dummy" rater equal to 1 for all +matches. Included so that \code{Euro2016} can be used directly with +\code{MCMCpack}'s models for pairwise comparisons.} + \item{team1}{The home team} + \item{team2}{The away team } + \item{winner}{The winner of the match. \code{NA} if a draw.} +} +} +\source{ +\url{https://en.wikipedia.org/wiki/UEFA_Euro_2016} +} +\description{ +Data on head-to-head outcomes from the 2016 UEFA European Football +Championship. +} +\keyword{datasets} diff -Nru r-cran-mcmcpack-1.5-0/man/HDPHMMnegbin.Rd r-cran-mcmcpack-1.6-0/man/HDPHMMnegbin.Rd --- r-cran-mcmcpack-1.5-0/man/HDPHMMnegbin.Rd 2021-01-19 05:26:43.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/HDPHMMnegbin.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -217,7 +217,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/HDPHMMpoisson.Rd r-cran-mcmcpack-1.6-0/man/HDPHMMpoisson.Rd --- r-cran-mcmcpack-1.5-0/man/HDPHMMpoisson.Rd 2021-01-19 05:26:43.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/HDPHMMpoisson.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -181,7 +181,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/HDPHSMMnegbin.Rd r-cran-mcmcpack-1.6-0/man/HDPHSMMnegbin.Rd --- r-cran-mcmcpack-1.5-0/man/HDPHSMMnegbin.Rd 2021-01-19 05:26:43.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/HDPHSMMnegbin.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -234,7 +234,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/HMMpanelFE.Rd r-cran-mcmcpack-1.6-0/man/HMMpanelFE.Rd --- r-cran-mcmcpack-1.5-0/man/HMMpanelFE.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/HMMpanelFE.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -213,5 +213,9 @@ Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.'' \emph{Journal of Econometrics}. 86: 221-241. + + Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. + ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical + Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. } \keyword{models} diff -Nru r-cran-mcmcpack-1.5-0/man/HMMpanelRE.Rd r-cran-mcmcpack-1.6-0/man/HMMpanelRE.Rd --- r-cran-mcmcpack-1.5-0/man/HMMpanelRE.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/HMMpanelRE.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -264,5 +264,9 @@ Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.'' \emph{Journal of Econometrics}. 86: 221-241. + +Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. +``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of +Statistical Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. } \keyword{models} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCbinaryChange.Rd r-cran-mcmcpack-1.6-0/man/MCMCbinaryChange.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCbinaryChange.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCbinaryChange.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -149,7 +149,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. 42(9): 1-21. -\url{https://www.jstatsoft.org/v42/i09/}. +\doi{10.18637/jss.v042.i09}. Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' \emph{Journal of the American Statistical diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCdynamicEI.Rd r-cran-mcmcpack-1.6-0/man/MCMCdynamicEI.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCdynamicEI.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCdynamicEI.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -206,7 +206,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Radford Neal. 2003. ``Slice Sampling" (with discussion). \emph{Annals of Statistics}, 31: 705-767. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCdynamicIRT1d.Rd r-cran-mcmcpack-1.6-0/man/MCMCdynamicIRT1d.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCdynamicIRT1d.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCdynamicIRT1d.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -302,7 +302,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. } \seealso{ \code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCfactanal.Rd r-cran-mcmcpack-1.6-0/man/MCMCfactanal.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCfactanal.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCfactanal.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -200,7 +200,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMChierEI.Rd r-cran-mcmcpack-1.6-0/man/MCMChierEI.Rd --- r-cran-mcmcpack-1.5-0/man/MCMChierEI.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMChierEI.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -174,7 +174,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCirt1d.Rd r-cran-mcmcpack-1.6-0/man/MCMCirt1d.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCirt1d.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCirt1d.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -213,7 +213,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCirtHier1d.Rd r-cran-mcmcpack-1.6-0/man/MCMCirtHier1d.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCirtHier1d.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCirtHier1d.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -288,7 +288,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCirtKd.Rd r-cran-mcmcpack-1.6-0/man/MCMCirtKd.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCirtKd.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCirtKd.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -215,7 +215,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCirtKdRob.Rd r-cran-mcmcpack-1.6-0/man/MCMCirtKdRob.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCirtKdRob.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCirtKdRob.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -335,7 +335,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMClogit.Rd r-cran-mcmcpack-1.6-0/man/MCMClogit.Rd --- r-cran-mcmcpack-1.5-0/man/MCMClogit.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMClogit.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -188,7 +188,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCmetrop1R.Rd r-cran-mcmcpack-1.6-0/man/MCMCmetrop1R.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCmetrop1R.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCmetrop1R.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -220,7 +220,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCmixfactanal.Rd r-cran-mcmcpack-1.6-0/man/MCMCmixfactanal.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCmixfactanal.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCmixfactanal.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -270,7 +270,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCmnl.Rd r-cran-mcmcpack-1.6-0/man/MCMCmnl.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCmnl.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCmnl.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -255,7 +255,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCnegbinChange.Rd r-cran-mcmcpack-1.6-0/man/MCMCnegbinChange.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCnegbinChange.Rd 2021-01-19 05:26:43.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCnegbinChange.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -218,7 +218,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCnegbin.Rd r-cran-mcmcpack-1.6-0/man/MCMCnegbin.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCnegbin.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCnegbin.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -160,7 +160,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCoprobitChange.Rd r-cran-mcmcpack-1.6-0/man/MCMCoprobitChange.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCoprobitChange.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCoprobitChange.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -194,7 +194,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.'' \emph{Journal of Econometrics}. 86: 221-241. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCoprobit.Rd r-cran-mcmcpack-1.6-0/man/MCMCoprobit.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCoprobit.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCoprobit.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -158,7 +158,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}. Springer: New York. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCordfactanal.Rd r-cran-mcmcpack-1.6-0/man/MCMCordfactanal.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCordfactanal.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCordfactanal.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -200,7 +200,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCpaircompare2dDP.Rd r-cran-mcmcpack-1.6-0/man/MCMCpaircompare2dDP.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCpaircompare2dDP.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCpaircompare2dDP.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -0,0 +1,351 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MCMCpaircompare2dDP.R +\name{MCMCpaircompare2dDP} +\alias{MCMCpaircompare2dDP} +\title{Markov Chain Monte Carlo for the Two-Dimensional Pairwise Comparisons +Model with Dirichlet Process Prior in Yu and Quinn (2021)} +\usage{ +MCMCpaircompare2dDP( + pwc.data, + theta.constraints = list(), + burnin = 1000, + mcmc = 20000, + thin = 1, + verbose = 0, + seed = NA, + gamma.start = NA, + theta.start = NA, + store.theta = TRUE, + store.gamma = FALSE, + tune = 0.3, + procrustes = FALSE, + alpha.start = 1, + cluster.max = 100, + cluster.mcmc = 500, + alpha.fixed = TRUE, + a0 = 1, + b0 = 1, + ... +) +} +\arguments{ +\item{pwc.data}{A data.frame containing the pairwise comparisons data. +Each row of \code{pwc.data} corresponds to a single pairwise comparison. +\code{pwc.data} needs to have exactly four columns. The first column +contains a unique identifier for the rater. Column two contains the unique +identifier for the first item being compared. Column three contains the +unique identifier for the second item being compared. Column four contains +the unique identifier of the item selected from the two items being +compared. If a tie occurred, the entry in the fourth column should be NA. +\strong{The identifiers in columns 2 through 4 must start with a letter. +Examples are provided below.}} + +\item{theta.constraints}{A list specifying possible simple equality or +inequality constraints on the item parameters. A + typical entry in the list has one of three forms: + \code{itemname=list(d,c)} which will constrain the dth dimension of + theta for the item named \code{itemname} to be equal to c, + \code{itemname=list(d,"+")} which will constrain the dth dimension of + theta for the item named \code{itemname} to be positive, and + \code{itemname=list(d, "-")} which will constrain the dth dimension of + theta for the item named \code{itemname} to be negative.} + +\item{burnin}{The number of burn-in iterations for the sampler.} + +\item{mcmc}{The number of Gibbs iterations for the sampler.} + +\item{thin}{The thinning interval used in the simulation. The number of +Gibbs iterations must be divisible by this value.} + +\item{verbose}{A switch which determines whether or not the progress of the +sampler is printed to the screen. If \code{verbose} is greater than 0 +output is printed to the screen every +\code{verbose}th iteration.} + +\item{seed}{The seed for the random number generator. If NA, the Mersenne +Twister generator is used with default seed 12345; if an integer is passed +it is used to seed the Mersenne twister. The user can also pass a list of +length two to use the L'Ecuyer random number generator, which is suitable +for parallel computation. The first element of the list is the L'Ecuyer +seed, which is a vector of length six or NA (if NA a default seed of +\code{rep(12345,6)} is used). The second element of list is a positive +substream number. See the MCMCpack specification for more details.} + +\item{gamma.start}{The starting value for the gamma vector. This +can either be a scalar or a column vector with dimension equal to the number +of raters. If this takes a scalar value, then that value will serve as the +starting value for all of the gammas. The default value of NA will set the +starting value of each gamma parameter to \eqn{\pi/4}.} + +\item{theta.start}{Starting values for the theta. Can be either a numeric +scalar, a J by 2 matrix (where J is the number of items compared), or NA. +If a scalar, all theta values are set to that value (except elements already +specified via theta.contraints. If NA, then non constrained elements of +theta are set equal to 0, elements constrained to be positive are set equal +to 0.5, elements constrained to be negative are set equal to -0.5 and +elements with equality constraints are set to satisfy those constraints.} + +\item{store.theta}{Should the theta draws be returned? Default is TRUE.} + +\item{store.gamma}{Should the gamma draws be returned? Default is TRUE.} + +\item{tune}{Tuning parameter for the random walk Metropolis proposal for +each gamma_i. \code{tune} is the width of the uniform proposal centered at +the current value of gamma_i. Must be a positive scalar.} + +\item{procrustes}{Should the theta and gamma draws be post-processed with +a Procrustes transformation? Default is FALSE. The Procrustes target matrix +is derived from the constrained elements of theta. Each row of theta that +has both theta values constrained is part of the of the target matrix. +Elements with equality constraints are set to those values. Elements +constrained to be positive are set to 1. Elements constrained to be negative +are set to -1. If \code{procrustes} is set to \code{TRUE} theta.constraints +must be set so that there are at least three rows of theta that have both +elements of theta constrained.} + +\item{alpha.start}{The starting value for the DP concentration parameter +alpha. Must be a positive scalar. Defaults to 1. If \code{alpha.fixed} is +set equal to TRUE, then alpha is held fixed at alpha.start.} + +\item{cluster.max}{The maximum number of clusters allowed in the +approximation to the DP prior for gamma. Defaults to 100. Must be a +positive integer.} + +\item{cluster.mcmc}{The number of additional MCMC iterations that are done +to sample each cluster-specific gamma value within one main MCMC iteration. +Must be a positive integer. Defaults to 500. Setting this to a lower value +speeds runtime at the cost of (possibly) worse mixing.} + +\item{alpha.fixed}{Logical value indicating whether the DP concentration +parameter alpha be held fixed (\code{TRUE}) or estimated (\code{FALSE}).} + +\item{a0}{The shape parameter of the gamma prior for alpha. This is the +same parameterization of the gamma distribution as R's internal +\code{rgamma()} function. Only relevant if \code{alpha.fixed} is set equal +to \code{FALSE}. Defaults to 1.} + +\item{b0}{The rate parameter of the gamma prior for alpha. This is the +same parameterization of the gamma distribution as R's internal +\code{rgamma()} function. Only relevant if \code{alpha.fixed} is set equal +to \code{FALSE}. Defaults to 1.} + +\item{...}{further arguments to be passed} +} +\value{ +An mcmc object that contains the posterior sample. This object can +be summarized by functions provided by the coda package. Most of the column +names of the mcmc object are self explanatory. Note however that the columns +with names of the form "cluster.[raterID]" give the cluster membership of +each rater at each stored MCMC iteration. Because of the possibility of +label switching, the particular values of these cluster membership variables +are not meaningful. What is meaningful is whether two raters share the same +cluster membership value at a particular MCMC iteration. This indicates +that those two raters were clustered together during that iteration. +Finally, note that the "n.clusters" column gives the number of distinct +gamma values at each iteration, i.e. the number of clusters at that +iteration. +} +\description{ +This function generates a sample from the posterior distribution of a +model for pairwise comparisons data with a probit link. Unlike standard +models for pairwise comparisons data, in this model the latent attribute +of each item being compared is a vector in two-dimensional Euclidean space. +} +\details{ +\code{MCMCpaircompare2d} uses the data augmentation approach of Albert and +Chib (1993) in conjunction with Gibbs and Metropolis-within-Gibbs steps +to fit the model. The user supplies data and a sample from the +posterior is returned as an \code{mcmc} object, which can be subsequently +analyzed in the \code{coda} package. + +The simulation is done in compiled C++ code to maximize efficiency. + +Please consult the \code{coda} package documentation for a comprehensive +list of functions that can be used to analyze the posterior sample. + +The model takes the following form: + +\deqn{i = 1,...,I \ \ \ \ (raters) } +\deqn{j = 1,...,J \ \ \ \ (items) } +\deqn{Y_{ijj'} = 1 \ \ if \ \ i \ \ chooses \ \ j \ \ over \ \ j'} +\deqn{Y_{ijj'} = 0 \ \ if \ \ i \ \ chooses \ \ j' \ \ over \ \ j} +\deqn{Y_{ijj'} = NA \ \ if \ \ i \ \ chooses \ \ neither} + +\deqn{\Pr(Y_{ijj'} = 1) = \Phi( \mathbf{z}_{i}' [\boldsymbol{\theta}_{j} - +\boldsymbol{\theta}_{ j'} ])} +\deqn{\mathbf{z}_{i}=[\cos(\gamma_{i}), \ \sin(\gamma_{i})]' } + +The following priors are assumed: +\deqn{\gamma_i \sim G} +\deqn{G \sim \mathcal{DP}(\alpha G_0)} +\deqn{G_0 = \mathcal{U}nif(0, \pi/2)} +\deqn{\alpha \sim \mathcal{G}amma(a_0, b_0)} +\deqn{\boldsymbol{\theta}_j \sim +\mathcal{N}_{2}(\mathbf{0}, \mathbf{I}_{2})} +For identification, some \eqn{\boldsymbol{\theta}_j}s are truncated +above or below 0, or fixed to constants. +} +\examples{ +\dontrun{ +## a synthetic data example +set.seed(123) + +I <- 65 ## number of raters +J <- 50 ## number of items to be compared + +## 3 clusters: +## raters 1 to 5 put most weight on dimension 1 +## raters 6 to 10 put most weight on dimension 2 +## raters 11 to I put substantial weight on both dimensions +gamma.true <- c(rep(0.05, 5), + rep(1.50, 5), + rep(0.7, I-10) ) +theta1.true <- rnorm(J, m=0, s=1) +theta2.true <- rnorm(J, m=0, s=1) +theta1.true[1] <- 2 +theta2.true[1] <- 2 +theta1.true[2] <- -2 +theta2.true[2] <- -2 +theta1.true[3] <- 2 +theta2.true[3] <- -2 + + + +n.comparisons <- 125 ## number of pairwise comparisons for each rater + +## generate synthetic data according to the assumed model +rater.id <- NULL +item.1.id <- NULL +item.2.id <- NULL +choice.id <- NULL +for (i in 1:I){ + for (c in 1:n.comparisons){ + rater.id <- c(rater.id, i+100) + item.numbers <- sample(1:J, size=2, replace=FALSE) + item.1 <- item.numbers[1] + item.2 <- item.numbers[2] + item.1.id <- c(item.1.id, item.1) + item.2.id <- c(item.2.id, item.2) + z <- c(cos(gamma.true[i]), sin(gamma.true[i])) + eta <- z[1] * (theta1.true[item.1] - theta1.true[item.2]) + + z[2] * (theta2.true[item.1] - theta2.true[item.2]) + prob.item.1.chosen <- pnorm(eta) + u <- runif(1) + if (u <= prob.item.1.chosen){ + choice.id <- c(choice.id, item.1) + } + else{ + choice.id <- c(choice.id, item.2) + } + } +} +item.1.id <- paste("item", item.1.id+100, sep=".") +item.2.id <- paste("item", item.2.id+100, sep=".") +choice.id <- paste("item", choice.id+100, sep=".") + +sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id) + + +## fit the model (should be run for more than 10500 iterations) +posterior <- MCMCpaircompare2dDP(pwc.data=sim.data, + theta.constraints=list(item.101=list(1,2), + item.101=list(2,2), + item.102=list(1,-2), + item.102=list(2,-2), + item.103=list(1,"+"), + item.103=list(2,"-")), + verbose=100, + burnin=500, mcmc=10000, thin=5, + cluster.mcmc=10, + store.theta=TRUE, store.gamma=TRUE, + tune=0.1) + + + + + +theta1.draws <- posterior[, grep("theta1", colnames(posterior))] +theta2.draws <- posterior[, grep("theta2", colnames(posterior))] +gamma.draws <- posterior[, grep("gamma", colnames(posterior))] + +theta1.post.med <- apply(theta1.draws, 2, median) +theta2.post.med <- apply(theta2.draws, 2, median) +gamma.post.med <- apply(gamma.draws, 2, median) + +theta1.post.025 <- apply(theta1.draws, 2, quantile, prob=0.025) +theta1.post.975 <- apply(theta1.draws, 2, quantile, prob=0.975) +theta2.post.025 <- apply(theta2.draws, 2, quantile, prob=0.025) +theta2.post.975 <- apply(theta2.draws, 2, quantile, prob=0.975) +gamma.post.025 <- apply(gamma.draws, 2, quantile, prob=0.025) +gamma.post.975 <- apply(gamma.draws, 2, quantile, prob=0.975) + + + +## compare estimates to truth +par(mfrow=c(2,2)) +plot(theta1.true, theta1.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +segments(x0=theta1.true, x1=theta1.true, + y0=theta1.post.025, y1=theta1.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + +plot(theta2.true, theta2.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +segments(x0=theta2.true, x1=theta2.true, + y0=theta2.post.025, y1=theta2.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + +plot(gamma.true, gamma.post.med, xlim=c(0, 1.6), ylim=c(0, 1.6), + col=rgb(0,0,0,0.3)) +segments(x0=gamma.true, x1=gamma.true, + y0=gamma.post.025, y1=gamma.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + + +## plot point estimates +plot(theta1.post.med, theta2.post.med, + xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +for (i in 1:length(gamma.post.med)){ + arrows(x0=0, y0=0, + x1=cos(gamma.post.med[i]), + y1=sin(gamma.post.med[i]), + col=rgb(1,0,0,0.2), len=0.05, lwd=0.5) +} + +} +} +\references{ +Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary +and Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, +669-679 + +Yu, Qiushi and Kevin M. Quinn. 2021. ``A Multidimensional Pairwise +Comparison Model for Heterogeneous Perceptions with an Application to +Modeling the Perceived Truthfulness of Public Statements on COVID-19.'' +University of Michigan Working Paper. + +Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: +Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. + +Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe +Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. + +Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output +Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. +\url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. +} +\seealso{ +\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, +\code{\link[MCMCpack]{MCMCpaircompare}}, +\code{\link[MCMCpack]{MCMCpaircompare2dDP}} +} +\author{ +Qiushi Yu and +Kevin M. Quinn +} +\keyword{models} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCpaircompare2d.Rd r-cran-mcmcpack-1.6-0/man/MCMCpaircompare2d.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCpaircompare2d.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCpaircompare2d.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -0,0 +1,304 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MCMCpaircompare2d.R +\name{MCMCpaircompare2d} +\alias{MCMCpaircompare2d} +\title{Markov Chain Monte Carlo for the Two-Dimensional Pairwise Comparisons +Model in Yu and Quinn (2021)} +\usage{ +MCMCpaircompare2d( + pwc.data, + theta.constraints = list(), + burnin = 1000, + mcmc = 20000, + thin = 1, + verbose = 0, + seed = NA, + gamma.start = NA, + theta.start = NA, + store.theta = TRUE, + store.gamma = TRUE, + tune = 0.3, + procrustes = FALSE, + ... +) +} +\arguments{ +\item{pwc.data}{A data.frame containing the pairwise comparisons data. +Each row of \code{pwc.data} corresponds to a single pairwise comparison. +\code{pwc.data} needs to have exactly four columns. The first column +contains a unique identifier for the rater. Column two contains the unique +identifier for the first item being compared. Column three contains the +unique identifier for the second item being compared. Column four contains +the unique identifier of the item selected from the two items being +compared. If a tie occurred, the entry in the fourth column should be NA. +\strong{The identifiers in columns 2 through 4 must start with a letter. +Examples are provided below.}} + +\item{theta.constraints}{A list specifying possible simple equality or +inequality constraints on the item parameters. A + typical entry in the list has one of three forms: + \code{itemname=list(d,c)} which will constrain the dth dimension of + theta for the item named \code{itemname} to be equal to c, + \code{itemname=list(d,"+")} which will constrain the dth dimension of + theta for the item named \code{itemname} to be positive, and + \code{itemname=list(d, "-")} which will constrain the dth dimension of + theta for the item named \code{itemname} to be negative.} + +\item{burnin}{The number of burn-in iterations for the sampler.} + +\item{mcmc}{The number of Gibbs iterations for the sampler.} + +\item{thin}{The thinning interval used in the simulation. The number of +Gibbs iterations must be divisible by this value.} + +\item{verbose}{A switch which determines whether or not the progress of the +sampler is printed to the screen. If \code{verbose} is greater than 0 +output is printed to the screen every +\code{verbose}th iteration.} + +\item{seed}{The seed for the random number generator. If NA, the Mersenne +Twister generator is used with default seed 12345; if an integer is passed +it is used to seed the Mersenne twister. The user can also pass a list of +length two to use the L'Ecuyer random number generator, which is suitable +for parallel computation. The first element of the list is the L'Ecuyer +seed, which is a vector of length six or NA (if NA a default seed of +\code{rep(12345,6)} is used). The second element of list is a positive +substream number. See the MCMCpack specification for more details.} + +\item{gamma.start}{The starting value for the gamma vector. This +can either be a scalar or a column vector with dimension equal to the number +of raters. If this takes a scalar value, then that value will serve as the +starting value for all of the gammas. The default value of NA will set the +starting value of each gamma parameter to \eqn{\pi/4}.} + +\item{theta.start}{Starting values for the theta. Can be either a numeric +scalar, a J by 2 matrix (where J is the number of items compared), or NA. +If a scalar, all theta values are set to that value (except elements already +specified via theta.contraints. If NA, then non constrained elements of +theta are set equal to 0, elements constrained to be positive are set equal +to 0.5, elements constrained to be negative are set equal to -0.5 and +elements with equality constraints are set to satisfy those constraints.} + +\item{store.theta}{Should the theta draws be returned? Default is TRUE.} + +\item{store.gamma}{Should the gamma draws be returned? Default is TRUE.} + +\item{tune}{Tuning parameter for the random walk Metropolis proposal for +each gamma_i. \code{tune} is the width of the uniform proposal centered at +the current value of gamma_i. Must be a positive scalar.} + +\item{procrustes}{Should the theta and gamma draws be post-processed with +a Procrustes transformation? Default is FALSE. The Procrustes target matrix +is derived from the constrained elements of theta. Each row of theta that +has both theta values constrained is part of the of the target matrix. +Elements with equality constraints are set to those values. Elements +constrained to be positive are set to 1. Elements constrained to be negative +are set to -1. If \code{procrustes} is set to \code{TRUE} theta.constraints +must be set so that there are at least three rows of theta that have both +elements of theta constrained.} + +\item{...}{further arguments to be passed} +} +\value{ +An mcmc object that contains the posterior sample. This object can +be summarized by functions provided by the coda package. +} +\description{ +This function generates a sample from the posterior distribution of a +model for pairwise comparisons data with a probit link. Unlike standard +models for pairwise comparisons data, in this model the latent attribute +of each item being compared is a vector in two-dimensional Euclidean space. +} +\details{ +\code{MCMCpaircompare2d} uses the data augmentation approach of Albert and +Chib (1993) in conjunction with Gibbs and Metropolis-within-Gibbs steps +to fit the model. The user supplies data and a sample from the +posterior is returned as an \code{mcmc} object, which can be subsequently +analyzed in the \code{coda} package. + +The simulation is done in compiled C++ code to maximize efficiency. + +Please consult the \code{coda} package documentation for a comprehensive +list of functions that can be used to analyze the posterior sample. + +The model takes the following form: + +\deqn{i = 1,...,I \ \ \ \ (raters) } +\deqn{j = 1,...,J \ \ \ \ (items) } +\deqn{Y_{ijj'} = 1 \ \ if \ \ i \ \ chooses \ \ j \ \ over \ \ j'} +\deqn{Y_{ijj'} = 0 \ \ if \ \ i \ \ chooses \ \ j' \ \ over \ \ j} +\deqn{Y_{ijj'} = NA \ \ if \ \ i \ \ chooses \ \ neither} + +\deqn{\Pr(Y_{ijj'} = 1) = \Phi( \mathbf{z}_{i}' [\boldsymbol{\theta}_{j} - +\boldsymbol{\theta}_{ j'} ])} +\deqn{\mathbf{z}_{i}=[\cos(\gamma_{i}), \ \sin(\gamma_{i})]' } + +The following priors are assumed: +\deqn{\gamma_i \sim \mathcal{U}nif(0, \ \pi/2)} +\deqn{\boldsymbol{\theta}_j \sim +\mathcal{N}_{2}(\mathbf{0}, \mathbf{I}_{2})} +For identification, some \eqn{\boldsymbol{\theta}_j}s are truncated +above or below 0, or fixed to constants. +} +\examples{ + + \dontrun{ +## a synthetic data example +set.seed(123) + +I <- 65 ## number of raters +J <- 50 ## number of items to be compared + + +## raters 1 to 5 put most weight on dimension 1 +## raters 6 to 10 put most weight on dimension 2 +## raters 11 to I put substantial weight on both dimensions +gamma.true <- c(runif(5, 0, 0.1), + runif(5, 1.47, 1.57), + runif(I-10, 0.58, 0.98) ) +theta1.true <- rnorm(J, m=0, s=1) +theta2.true <- rnorm(J, m=0, s=1) +theta1.true[1] <- 2 +theta2.true[1] <- 2 +theta1.true[2] <- -2 +theta2.true[2] <- -2 +theta1.true[3] <- 2 +theta2.true[3] <- -2 + + + +n.comparisons <- 125 ## number of pairwise comparisons for each rater + +## generate synthetic data according to the assumed model +rater.id <- NULL +item.1.id <- NULL +item.2.id <- NULL +choice.id <- NULL +for (i in 1:I){ + for (c in 1:n.comparisons){ + rater.id <- c(rater.id, i+100) + item.numbers <- sample(1:J, size=2, replace=FALSE) + item.1 <- item.numbers[1] + item.2 <- item.numbers[2] + item.1.id <- c(item.1.id, item.1) + item.2.id <- c(item.2.id, item.2) + z <- c(cos(gamma.true[i]), sin(gamma.true[i])) + eta <- z[1] * (theta1.true[item.1] - theta1.true[item.2]) + + z[2] * (theta2.true[item.1] - theta2.true[item.2]) + prob.item.1.chosen <- pnorm(eta) + u <- runif(1) + if (u <= prob.item.1.chosen){ + choice.id <- c(choice.id, item.1) + } + else{ + choice.id <- c(choice.id, item.2) + } + } +} +item.1.id <- paste("item", item.1.id+100, sep=".") +item.2.id <- paste("item", item.2.id+100, sep=".") +choice.id <- paste("item", choice.id+100, sep=".") + +sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id) + + +## fit the model +posterior <- MCMCpaircompare2d(pwc.data=sim.data, + theta.constraints=list(item.101=list(1,2), + item.101=list(2,2), + item.102=list(1,-2), + item.102=list(2,-2), + item.103=list(1,"+"), + item.103=list(2,"-")), + verbose=1000, + burnin=500, mcmc=20000, thin=10, + store.theta=TRUE, store.gamma=TRUE, tune=0.5) + + + + + +theta1.draws <- posterior[, grep("theta1", colnames(posterior))] +theta2.draws <- posterior[, grep("theta2", colnames(posterior))] +gamma.draws <- posterior[, grep("gamma", colnames(posterior))] + +theta1.post.med <- apply(theta1.draws, 2, median) +theta2.post.med <- apply(theta2.draws, 2, median) +gamma.post.med <- apply(gamma.draws, 2, median) + +theta1.post.025 <- apply(theta1.draws, 2, quantile, prob=0.025) +theta1.post.975 <- apply(theta1.draws, 2, quantile, prob=0.975) +theta2.post.025 <- apply(theta2.draws, 2, quantile, prob=0.025) +theta2.post.975 <- apply(theta2.draws, 2, quantile, prob=0.975) +gamma.post.025 <- apply(gamma.draws, 2, quantile, prob=0.025) +gamma.post.975 <- apply(gamma.draws, 2, quantile, prob=0.975) + + + +## compare estimates to truth +par(mfrow=c(2,2)) +plot(theta1.true, theta1.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +segments(x0=theta1.true, x1=theta1.true, + y0=theta1.post.025, y1=theta1.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + +plot(theta2.true, theta2.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +segments(x0=theta2.true, x1=theta2.true, + y0=theta2.post.025, y1=theta2.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + +plot(gamma.true, gamma.post.med, xlim=c(0, 1.6), ylim=c(0, 1.6), + col=rgb(0,0,0,0.3)) +segments(x0=gamma.true, x1=gamma.true, + y0=gamma.post.025, y1=gamma.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + + +## plot point estimates +plot(theta1.post.med, theta2.post.med, + xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +for (i in 1:length(gamma.post.med)){ + arrows(x0=0, y0=0, + x1=cos(gamma.post.med[i]), + y1=sin(gamma.post.med[i]), + col=rgb(1,0,0,0.2), len=0.05, lwd=0.5) +} +} +} +\references{ +Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary +and Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, +669-679 + +Yu, Qiushi and Kevin M. Quinn. 2021. ``A Multidimensional Pairwise +Comparison Model for Heterogeneous Perceptions with an Application to +Modeling the Perceived Truthfulness of Public Statements on COVID-19.'' +University of Michigan Working Paper. + +Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: +Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. + +Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe +Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. + +Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output +Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. +\url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. +} +\seealso{ +\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, +\code{\link[MCMCpack]{MCMCpaircompare}}, +\code{\link[MCMCpack]{MCMCpaircompare2dDP}} +} +\author{ +Qiushi Yu and +Kevin M. Quinn +} +\keyword{models} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCpaircompare.Rd r-cran-mcmcpack-1.6-0/man/MCMCpaircompare.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCpaircompare.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCpaircompare.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -0,0 +1,266 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/MCMCpaircompare.R +\name{MCMCpaircompare} +\alias{MCMCpaircompare} +\title{Markov Chain Monte Carlo for a Pairwise Comparisons Model with Probit Link} +\usage{ +MCMCpaircompare( + pwc.data, + theta.constraints = list(), + alpha.fixed = FALSE, + burnin = 1000, + mcmc = 20000, + thin = 1, + verbose = 0, + seed = NA, + alpha.start = NA, + a = 0, + A = 0.25, + store.theta = TRUE, + store.alpha = FALSE, + ... +) +} +\arguments{ +\item{pwc.data}{A data.frame containing the pairwise comparisons data. +Each row of \code{pwc.data} corresponds to a single pairwise comparison. +\code{pwc.data} needs to have exactly four columns. The first column +contains a unique identifier for the rater. Column two contains the unique +identifier for the first item being compared. Column three contains the +unique identifier for the second item being compared. Column four contains +the unique identifier of the item selected from the two items being +compared. If a tie occurred, the entry in the fourth column should be NA. +For applications without raters (such as sports competitions) all entries +in the first column should be set to a single value and \code{alpha.fixed} +(see below) should be set to \code{TRUE}. \strong{The identifiers in +columns 2 through 4 must start with a letter. Examples are provided below.}} + +\item{theta.constraints}{A list specifying possible simple equality or +inequality constraints on the item parameters. A typical entry in the +list has one of three forms: \code{itemname=c} which will constrain the +item parameter for the item named \code{itemname} to be equal to c, +\code{itemname="+"} which will constrain the item parameter for the +item named \code{itemname} to be positive, and \code{itemname="-"} which +will constrain the item parameter for the item named \code{itemname} to +be negative.} + +\item{alpha.fixed}{Should alpha be fixed to a constant value of 1 for all +raters? Default is FALSE. If set to FALSE, an alpha value is estimated for +each rater.} + +\item{burnin}{The number of burn-in iterations for the sampler.} + +\item{mcmc}{The number of Gibbs iterations for the sampler.} + +\item{thin}{The thinning interval used in the simulation. The number of +Gibbs iterations must be divisible by this value.} + +\item{verbose}{A switch which determines whether or not the progress of the +sampler is printed to the screen. If \code{verbose} is greater than 0 +output is printed to the screen every +\code{verbose}th iteration.} + +\item{seed}{The seed for the random number generator. If NA, the Mersenne +Twister generator is used with default seed 12345; if an integer is passed +it is used to seed the Mersenne twister. The user can also pass a list of +length two to use the L'Ecuyer random number generator, which is suitable +for parallel computation. The first element of the list is the L'Ecuyer +seed, which is a vector of length six or NA (if NA a default seed of +\code{rep(12345,6)} is used). The second element of list is a positive +substream number. See the MCMCpack specification for more details.} + +\item{alpha.start}{The starting value for the alpha vector. This +can either be a scalar or a column vector with dimension equal to the number +of alphas. If this takes a scalar value, then that value will serve as the +starting value for all of the alphas. The default value of NA will set the +starting value of each alpha parameter to 1.} + +\item{a}{The prior mean of alpha. Must be a scalar. Default is 0.} + +\item{A}{The prior precision of alpha. Must be a positive scalar. +Default is 0.25 (prior variance is 4).} + +\item{store.theta}{Should the theta draws be returned? Default is TRUE.} + +\item{store.alpha}{Should the alpha draws be returned? Default is FALSE.} + +\item{...}{further arguments to be passed} +} +\value{ +An mcmc object that contains the posterior sample. This object can +be summarized by functions provided by the coda package. +} +\description{ +This function generates a sample from the posterior distribution of a +model for pairwise comparisons data with a probit link. Thurstone's model +is a special case of this model when the \eqn{\alpha} parameter is fixed at +1. +} +\details{ +\code{MCMCpaircompare} uses the data augmentation approach of Albert and +Chib (1993). The user supplies data and priors, and a sample from the +posterior is returned as an \code{mcmc} object, which can be subsequently +analyzed in the \code{coda} package. + +The simulation is done in compiled C++ code to maximize efficiency. + +Please consult the \code{coda} package documentation for a comprehensive +list of functions that can be used to analyze the posterior sample. + +The model takes the following form: + +\deqn{i = 1,...,I \ \ \ \ (raters) } +\deqn{j = 1,...,J \ \ \ \ (items) } +\deqn{Y_{ijj'} = 1 \ \ if \ \ i \ \ chooses \ \ j \ \ over \ \ j'} +\deqn{Y_{ijj'} = 0 \ \ if \ \ i \ \ chooses \ \ j' \ \ over \ \ j} +\deqn{Y_{ijj'} = NA \ \ if \ \ i \ \ chooses \ \ neither} + +\deqn{Pr(Y_{ijj'} = 1) = \Phi( \alpha_{i} [\theta_{j} - \theta_{ j'} ] ) } + +The following Gaussian priors are assumed: +\deqn{\alpha_i \sim \mathcal{N}(a, A^{-1})} +\deqn{\theta_j \sim \mathcal{N}(0, 1)} +For identification, some \eqn{\theta_j}s are truncated above or below 0, +or fixed to constants. +} +\examples{ + + \dontrun{ + ## Euro 2016 example + data(Euro2016) + +posterior1 <- MCMCpaircompare(pwc.data=Euro2016, + theta.constraints=list(Ukraine="-", + Portugal="+"), + alpha.fixed=TRUE, + verbose=10000, + burnin=10000, mcmc=500000, thin=100, + store.theta=TRUE, store.alpha=FALSE) + +## alternative identification constraints +posterior2 <- MCMCpaircompare(pwc.data=Euro2016, + theta.constraints=list(Ukraine="-", + Portugal=1), + alpha.fixed=TRUE, + verbose=10000, + burnin=10000, mcmc=500000, thin=100, + store.theta=TRUE, store.alpha=FALSE) + + + + + + + + +## a synthetic data example with estimated rater-specific parameters +set.seed(123) + +I <- 65 ## number of raters +J <- 50 ## number of items to be compared + + +## raters 1 to 5 have less sensitivity to stimuli than raters 6 through I +alpha.true <- c(rnorm(5, m=0.2, s=0.05), rnorm(I - 5, m=1, s=0.1)) +theta.true <- sort(rnorm(J, m=0, s=1)) + +n.comparisons <- 125 ## number of pairwise comparisons for each rater + +## generate synthetic data according to the assumed model +rater.id <- NULL +item.1.id <- NULL +item.2.id <- NULL +choice.id <- NULL +for (i in 1:I){ + for (c in 1:n.comparisons){ + rater.id <- c(rater.id, i+100) + item.numbers <- sample(1:J, size=2, replace=FALSE) + item.1 <- item.numbers[1] + item.2 <- item.numbers[2] + item.1.id <- c(item.1.id, item.1) + item.2.id <- c(item.2.id, item.2) + eta <- alpha.true[i] * (theta.true[item.1] - theta.true[item.2]) + prob.item.1.chosen <- pnorm(eta) + u <- runif(1) + if (u <= prob.item.1.chosen){ + choice.id <- c(choice.id, item.1) + } + else{ + choice.id <- c(choice.id, item.2) + } + } +} +item.1.id <- paste("item", item.1.id+100, sep=".") +item.2.id <- paste("item", item.2.id+100, sep=".") +choice.id <- paste("item", choice.id+100, sep=".") + +sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id) + + +## fit the model +posterior <- MCMCpaircompare(pwc.data=sim.data, + theta.constraints=list(item.101=-2, + item.150=2), + alpha.fixed=FALSE, + verbose=10000, + a=0, A=0.5, + burnin=10000, mcmc=200000, thin=100, + store.theta=TRUE, store.alpha=TRUE) + +theta.draws <- posterior[, grep("theta", colnames(posterior))] +alpha.draws <- posterior[, grep("alpha", colnames(posterior))] + +theta.post.med <- apply(theta.draws, 2, median) +alpha.post.med <- apply(alpha.draws, 2, median) + +theta.post.025 <- apply(theta.draws, 2, quantile, prob=0.025) +theta.post.975 <- apply(theta.draws, 2, quantile, prob=0.975) +alpha.post.025 <- apply(alpha.draws, 2, quantile, prob=0.025) +alpha.post.975 <- apply(alpha.draws, 2, quantile, prob=0.975) + +## compare estimates to truth +par(mfrow=c(1,2)) +plot(theta.true, theta.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), + col=rgb(0,0,0,0.3)) +segments(x0=theta.true, x1=theta.true, + y0=theta.post.025, y1=theta.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + +plot(alpha.true, alpha.post.med, xlim=c(0, 1.2), ylim=c(0, 3), + col=rgb(0,0,0,0.3)) +segments(x0=alpha.true, x1=alpha.true, + y0=alpha.post.025, y1=alpha.post.975, + col=rgb(0,0,0,0.3)) +abline(0, 1, col=rgb(1,0,0,0.5)) + +} + +} +\references{ +Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary +and Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, +669-679 + +Yu, Qiushi and Kevin M. Quinn. 2021. ``A Multidimensional Pairwise +Comparison Model for Heterogeneous Perception with an Application to +Modeling the Perceived Truthfulness of Public Statements on COVID-19.'' +University of Michigan Working Paper. + +Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: +Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. + +Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe +Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. + +Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output +Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. +\url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. +} +\seealso{ +\code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, +\code{\link[MCMCpack]{MCMCpaircompare2d}}, +\code{\link[MCMCpack]{MCMCpaircompare2dDP}} +} +\keyword{models} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCpoissonChange.Rd r-cran-mcmcpack-1.6-0/man/MCMCpoissonChange.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCpoissonChange.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCpoissonChange.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -214,7 +214,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' \emph{Journal of the American Statistical Association}. 90: 1313-1321. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCpoisson.Rd r-cran-mcmcpack-1.6-0/man/MCMCpoisson.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCpoisson.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCpoisson.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -129,7 +129,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCprobitChange.Rd r-cran-mcmcpack-1.6-0/man/MCMCprobitChange.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCprobitChange.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCprobitChange.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -178,7 +178,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point models.'' \emph{Journal of Econometrics}. 86: 221-241. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCprobit.Rd r-cran-mcmcpack-1.6-0/man/MCMCprobit.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCprobit.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCprobit.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -136,7 +136,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +42(9): 1-21. \doi{10.18637/jss.v042.i09}. Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' \emph{Journal of the American Statistical Association}. 90: 1313-1321. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCregressChange.Rd r-cran-mcmcpack-1.6-0/man/MCMCregressChange.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCregressChange.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCregressChange.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -232,6 +232,10 @@ Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point models." \emph{Journal of Econometrics}. 86: 221-241. + +Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. +``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of +Statistical Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. } \seealso{ \code{\link{plotState}}, \code{\link{plotChangepoint}} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCregress.Rd r-cran-mcmcpack-1.6-0/man/MCMCregress.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCregress.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCregress.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -152,7 +152,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' \emph{Journal of the American Statistical Association}. 90: 1313-1321. diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCSVDreg.Rd r-cran-mcmcpack-1.6-0/man/MCMCSVDreg.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCSVDreg.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCSVDreg.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -148,7 +148,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. 42(9): 1-21. -\url{https://www.jstatsoft.org/v42/i09/}. +\doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} diff -Nru r-cran-mcmcpack-1.5-0/man/MCMCtobit.Rd r-cran-mcmcpack-1.6-0/man/MCMCtobit.Rd --- r-cran-mcmcpack-1.5-0/man/MCMCtobit.Rd 2021-01-19 04:08:54.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/man/MCMCtobit.Rd 2021-10-06 01:57:26.000000000 +0000 @@ -157,7 +157,7 @@ \references{ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/MD5 r-cran-mcmcpack-1.6-0/MD5 --- r-cran-mcmcpack-1.5-0/MD5 2021-01-20 11:50:12.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/MD5 2021-10-06 05:40:02.000000000 +0000 @@ -1,55 +1,58 @@ -658fbbabb8c88e15c2379df0c34475d9 *DESCRIPTION -4d49b8c5f1f90650c8967f18a5d1b2cc *NAMESPACE +12828356ea09f10e8788158d0b427708 *DESCRIPTION +5cdd3af9ac0f101a82634541b96c65c0 *NAMESPACE 489c13dcd0f8b0bf0109cb26e339c489 *R/BayesFactors.R -3b8376853c83eeb3891e17a928039ca5 *R/HDPHMMnegbin.R -97b87a802e67caf114c644ff38c6a2c5 *R/HDPHMMpoisson.R -cfcf448563212d5a0359bb60d7f1c4ef *R/HDPHSMMnegbin.R -a59379c34c30105026e3bb07bec203e9 *R/HMMpanelFE.R -130d29a6e4d5338965e4a84b508a4af8 *R/HMMpanelRE.R -74bb0fe4f190b8621d22d38ed99b7f36 *R/MCMCSVDreg.R -399cc1ddfd56b6aa0e5d4b07468ed142 *R/MCMCbinaryChange.R -12c7517346c676830f6168d8f3114ea1 *R/MCMCdynamicEI.R +a73077aedf968174f4b546177f872332 *R/HDPHMMnegbin.R +09c65a1b7b643e698b939bf2262baec2 *R/HDPHMMpoisson.R +883a4527d857b56a22aa0b158720c178 *R/HDPHSMMnegbin.R +a0679289499facd20ab89b3a9347df3c *R/HMMpanelFE.R +9693b70572cdc96a8749a11c5250cfab *R/HMMpanelRE.R +9dcbb57b01995bef885eeef81b161219 *R/MCMCSVDreg.R +245ad43374868e25246d42089935147d *R/MCMCbinaryChange.R +03a2a080d31f3558030719d2e30c827c *R/MCMCdynamicEI.R 9ef25b73768c661af8ec2f2dd651391d *R/MCMCdynamicIRT1d-b.R -d37595c66704e8c70efb5b3cddf3f46c *R/MCMCdynamicIRT1d.R -2406e062966b231149f97679678610e8 *R/MCMCfactanal.R +6c854948a5cd6fe44806861a138f7136 *R/MCMCdynamicIRT1d.R +b587f167cac6ec0dd5b305563f344374 *R/MCMCfactanal.R dcd83d013877c39ace4ca33c4357c779 *R/MCMChierBetaBinom.R -579085714aa3e3cd7b2d73e7923a9cd2 *R/MCMChierEI.R +da6125335720e4741ff2c27e47cbaf67 *R/MCMChierEI.R 84093773bf1337ed03ab323742acc01d *R/MCMChlogit.R ddaa0f30be43cda30f816a7dd07ce354 *R/MCMChpoisson.R dce402371cba577cd7967d40c78b7f17 *R/MCMChregress.R -3ad61744e514ae834f157b0915442970 *R/MCMCirt1d.R -1ee0d7e705dcff3024980d3d1bfba824 *R/MCMCirtHier1d.R -cba7eddfad83f73c3fb5752ed6807dec *R/MCMCirtKd.R -fb0e9d49cd287b83c7a155aa02f7fea9 *R/MCMCirtKdRob.R -7d54fc18718315f8ba9546b6e3057742 *R/MCMClogit.R -701a1be796657e418ccfb6a1c22a99bb *R/MCMCmetrop1R.R -6258c7cbd7c982e69cceff339cc3e0b9 *R/MCMCmixfactanal.R -fdbdedea5e54acdaba83e320ce89cab5 *R/MCMCmnl.R -096e2b86f97cfdf17bcdef5cc631b24d *R/MCMCnegbin.R -3a04b8f7f7830280882284675f5e963b *R/MCMCnegbinChange.R -3340379853bef5f3d7d535e342dc34fe *R/MCMCoprobit.R -80718bde46823b2081a4548898da8c0c *R/MCMCoprobitChange.R -fc127100e0a4c96e8447d73ce17d678f *R/MCMCordfactanal.R +83cb5431ecc4ff9ad5f5f3530e3b3484 *R/MCMCirt1d.R +a9cc3eaf15fee57b8d2d9188ce57a850 *R/MCMCirtHier1d.R +6ebe606753f4238e7f0a80bdff0fe96c *R/MCMCirtKd.R +c404a10ed9422cdde7eb3f04374330d9 *R/MCMCirtKdRob.R +9ff2828825a0b37f221e7394cfd83359 *R/MCMClogit.R +d94303bc04acd8dcbadb5fecc0bcaa1c *R/MCMCmetrop1R.R +35adb8c9f932bd128d5412d524ae1d9a *R/MCMCmixfactanal.R +62086c625aa86ddd8145cda2aedacc57 *R/MCMCmnl.R +6605b03005420d2954de4b184db08137 *R/MCMCnegbin.R +14eed04fbcf387b0f8c27a49048eb6f1 *R/MCMCnegbinChange.R +08a402d761d51a1ccb550ff5d43b8ac6 *R/MCMCoprobit.R +cd59ac31b578700118b41faf3562afb4 *R/MCMCoprobitChange.R +2f9a4e72a33836b981a7cc4bd65f8dad *R/MCMCordfactanal.R fb8501996a264375e1e9e6daeed449c3 *R/MCMCpack-package.R -7b4fa9b995673c70c228b441a9cb4722 *R/MCMCpoisson.R -7d599cfb26b407109c070d4a03fe01ad *R/MCMCpoissonChange.R -6a88755a655fee3149246ba7526705b8 *R/MCMCprobit.R -b75a93c40ef592b6d7b3a548534a3333 *R/MCMCprobitChange.R +dcbb50f17d47b3c53583921a2ecbf240 *R/MCMCpaircompare.R +71725c07d979f089e2a2d0070fd882a7 *R/MCMCpaircompare2d.R +9dc2d3a83782854142f9a1b4674f6170 *R/MCMCpaircompare2dDP.R +4071148d9ad71b51fa22a7031e74aa92 *R/MCMCpoisson.R +4b7879834b4788ba36acff4762003942 *R/MCMCpoissonChange.R +8894b049eac0195bb81a83199eef4f73 *R/MCMCprobit.R +dee7081e726e16ad62d66aa59871be69 *R/MCMCprobitChange.R e4b5208b09de1590187c700bd450bf12 *R/MCMCquantreg.R -9930c27c2cf3c57b132095fb21e7f7db *R/MCMCregress.R -04278767133395d08b4497ae8a59b757 *R/MCMCregressChange.R +59e6304502af102c3d2be21fcd5f8163 *R/MCMCregress.R +3d8cea6297ec9c1cd48bea3d162e384d *R/MCMCregressChange.R cfaa7b7b68eaa8ebb77f5964e376c5c1 *R/MCMCresidualBreakAnalysis.R -d5752d10b4451f10cdd23d3ad4ac0d2d *R/MCMCtobit.R +7b60264616e1c51efc61ccb3f78f35d4 *R/MCMCtobit.R 4afd92a950cb5580bb0d21f05332f3a6 *R/MCmodels.R cc61576fd01bc5f2c3bf3291ac85f24c *R/SSVSquantreg.R 860f50b8e3d79da3d6a28d95744fba41 *R/SSVSquantregsummary.R fdfcfe9909cadaf15c47bbb604be2257 *R/automate.R 80bedd199cd49c78355d3383c0170832 *R/btsutil.R -77bce14a19cc9cc19b0101ff87dbcdd2 *R/data.R +45caf8e800f202268a86d5c128129b1a *R/data.R d51c4bb909bbbde846dd62441303a9c1 *R/distn.R 1589cbe1b84ab5b397b96ef59ba79422 *R/hdp-utils.R 18daa9fee800c7d8d5827405738531e2 *R/hidden-hmodels.R -14032465f6fb89990e246248c10aa9d7 *R/hidden.R +bab89af4304df7fc21f018ac3fb3f0ca *R/hidden.R 9ac7fe5740a008b7d085f4e6702f0924 *R/make.breaklist.R d451566d6a37400a46248a807dffd0e8 *R/procrust.R 8b765f0882bbd8ba6ff3e19da454c4fc *R/scythe.R @@ -59,56 +62,62 @@ 175f5ed085a3f66d7cb87d64bac0e344 *R/utility.R 7e2021def8a88edb41fe78d0a0ae642c *R/zzz.R da48e08cde3dff063f68dbf54c82f29b *README +b711bc619e0ce0d41f279dfb504feb5a *build/partial.rdb 31f24d1e199bde63095dcfa4556e848b *cleanup ce8d6bc2b702996ba91799bc33b82840 *configure 21e1a006fb388ed39daf16e6d2ec9549 *configure.ac +bbfdf82e8e3c4ff7d755cd4077e3ed08 *data/Euro2016.rda e17202d22ce5ab5ab7b01e3b248a4008 *data/Nethvote.rda 79c0617133751c81adad3dbe643aafea *data/PErisk.rda 93e636ac58b966f9a74315f51be07d25 *data/Rehnquist.rda accb7b5b46207a02e49de2702a6faff4 *data/Senate.rda ce79b1b47b561cb44dd4f93589d3d755 *data/SupremeCourt.rda -c600a9bfe75dee1a0f39db92a1d34459 *inst/CITATION -75ca22638ae19756a9981bd0146931b0 *inst/HISTORY.txt +c94fffe49ba014cbc6708326fbda0364 *inst/CITATION +954d223935ee10145211f1f6ce0ae13f *inst/HISTORY.txt cb4856b66283c4e98efeee9b0866b4ce *man/BayesFactor.Rd 3fde04199db6c766bf4474f697236a44 *man/Dirichlet.Rd -83609172d12ce7cfd8d4468a57e121c1 *man/HDPHMMnegbin.Rd -dcf5e4f2f42c757c648b399f301cbfa8 *man/HDPHMMpoisson.Rd -8767b819d6ccca79cd0a009f357b4e02 *man/HDPHSMMnegbin.Rd -80c14fe410446a6c3b0d9f6f73fd3f02 *man/HMMpanelFE.Rd -14601223980bfabd06e385cb5161e213 *man/HMMpanelRE.Rd +558a86fab20e23525e33a6bb7e74ab88 *man/Euro2016.Rd +5c58e90243e4c1feb29bb083b606c999 *man/HDPHMMnegbin.Rd +7c3c52fdf3ea032ce2803901d6756cc6 *man/HDPHMMpoisson.Rd +e28ff313b8cc179071aec612a2cf5e45 *man/HDPHSMMnegbin.Rd +c85e91a0cb46c965edd5736cdc4be642 *man/HMMpanelFE.Rd +f5014fc4315922cc51c16a9866c87c1f *man/HMMpanelRE.Rd c662ff9b3f3383e49351e15face53c72 *man/InvGamma.Rd 6e36da2c16cc559dccc1edc194032a7d *man/InvWishart.Rd -24a1a42c8f383e8aaf5e824489de3b53 *man/MCMCSVDreg.Rd -03209378cfa20c47e1283139b4d8a6c2 *man/MCMCbinaryChange.Rd -873f960835584993be2373eff8abaffd *man/MCMCdynamicEI.Rd -023d090e5749eb9d520ab632837d3c94 *man/MCMCdynamicIRT1d.Rd -0bb8b6ad61b52179e464b345db25cb1e *man/MCMCfactanal.Rd -eb2fc576e01714bbb9d3a53fee18f639 *man/MCMChierEI.Rd +873a6919ce184d44df96517d663b0e02 *man/MCMCSVDreg.Rd +2908dba0b0b47d87ce3622ee61da9b76 *man/MCMCbinaryChange.Rd +b884c99098d20a0b921cdedfd388779a *man/MCMCdynamicEI.Rd +86010ae32c29de3b92ce9feeb17b235b *man/MCMCdynamicIRT1d.Rd +295545010c95efeee55d490a9b1cd356 *man/MCMCfactanal.Rd +da0c55a99f25ea4592fcfa913877250f *man/MCMChierEI.Rd cf552f93f0cda2161cc7e12a17a32445 *man/MCMChlogit.Rd 757a4181bf12c3c1a26d07ed7e4f2697 *man/MCMChpoisson.Rd 9f4b6af1947a221ef87b16b72b140804 *man/MCMChregress.Rd -5afe5a4d8ca28bc898e7d40626f1ae15 *man/MCMCirt1d.Rd -12cfa090eeff1cb00f15a42e94d6b2eb *man/MCMCirtHier1d.Rd -5131eafac454530a2b1db6927997f57f *man/MCMCirtKd.Rd -209ad56615a02b132f2273ff167fbcbf *man/MCMCirtKdRob.Rd -32b0fabae5632fd9729b05230ef6f359 *man/MCMClogit.Rd -8296d4388463e1213c6d34ebb322b7ec *man/MCMCmetrop1R.Rd -02b4a4a542c0291fa70dc58783a2c2d0 *man/MCMCmixfactanal.Rd -8b55ec42db6737dd0fbf8dc81086b03e *man/MCMCmnl.Rd -ed6b362cdd8060d799efc1e0b37d647d *man/MCMCnegbin.Rd -5c602f942c5d3bd92c8187404090f66b *man/MCMCnegbinChange.Rd -639498ea6b3308a770c730f507d3fd31 *man/MCMCoprobit.Rd -10a0df72f58d8167c3a0673a52c6c42a *man/MCMCoprobitChange.Rd -a1bc9c1787d370393a3d16e92e3768b2 *man/MCMCordfactanal.Rd -820f8f839e7cbb83542bda75949940ad *man/MCMCpoisson.Rd -bbc8ff0e66e544fbb97d7e96dea5e744 *man/MCMCpoissonChange.Rd -0f163ada278e4a9b0405080d0eea2e9e *man/MCMCprobit.Rd -eaac408840fa89dc3bd949b35b181d90 *man/MCMCprobitChange.Rd +937db4bc90138817186c9141b66cf194 *man/MCMCirt1d.Rd +2e39777b285157eb1f742b7d36c1e4dd *man/MCMCirtHier1d.Rd +62fe2b231bf5c57adddd824f1e12de93 *man/MCMCirtKd.Rd +448c5947a8fa5880b53b63ce70da2fed *man/MCMCirtKdRob.Rd +24ce6124e98a706b080e0336c453fdbb *man/MCMClogit.Rd +28daf408c1ad06a47a41cd081f262fc0 *man/MCMCmetrop1R.Rd +983ad537525d0ffa42d1f0ac8c31bb80 *man/MCMCmixfactanal.Rd +5d65e5d4ea252c89fba433d3d97765b1 *man/MCMCmnl.Rd +bd3b07be733f834230b509defee21b57 *man/MCMCnegbin.Rd +267bb85012a7eccd89174b6e105bda92 *man/MCMCnegbinChange.Rd +584d506ed31167dee991c5a59466e9f8 *man/MCMCoprobit.Rd +a46af578f6ddd82726bba9465ab93354 *man/MCMCoprobitChange.Rd +1d51000fb7ee1b713487005a6afbc7e1 *man/MCMCordfactanal.Rd +737450145c38d759d90b4b21d9d1ab81 *man/MCMCpaircompare.Rd +120aceadcf648c76eef4ebe6fa80585a *man/MCMCpaircompare2d.Rd +27eec23254770c5d3962fa4028596a7f *man/MCMCpaircompare2dDP.Rd +1a6124c049d8fb7d1d6b39a87d585251 *man/MCMCpoisson.Rd +be76b6e9f110598c3cf9632d6cbbe2dc *man/MCMCpoissonChange.Rd +67d906617a70b7ef5abc843f75274b38 *man/MCMCprobit.Rd +b07efc5aff5f3073624a263f796efa3b *man/MCMCprobitChange.Rd 4a9358259cfb9a1deeba3c0f840ea512 *man/MCMCquantreg.Rd -e699ac35a8384f2f07a0d2ffcc0e91a6 *man/MCMCregress.Rd -38e8e5660e213fd857fa56f47d33dd95 *man/MCMCregressChange.Rd +c12594b91f315501282c53cbb59ee5f9 *man/MCMCregress.Rd +4a024563a6c184c2b5d76da59c44eff2 *man/MCMCregressChange.Rd de9c9c0ab4050d1d1931b4d27a56093a *man/MCMCresidualBreakAnalysis.Rd -0cf7246353249ecd5ba5a3523d99dcd8 *man/MCMCtobit.Rd +59aa441e819ee5113e38fcc63ffbe45a *man/MCMCtobit.Rd da7568c4e553c121d5f953d941118076 *man/MCbinomialbeta.Rd 842c56f0dc25f01face9f7fb416dea53 *man/MCmultinomdirichlet.Rd d79e5091a65f40dab48903b9e19592e9 *man/MCnormalnormal.Rd @@ -142,7 +151,7 @@ 1d71348836f0c214c36bd87b0204a488 *man/xpnd.Rd 5ca10cd5c2050d7ae43a8d4756ab9006 *src/HMMmultivariateGaussian.cc 8dc624f6b2ea2967170da8cc98a04bdb *src/MCMCdynamicEI.cc -a35f29e283080796db85b4c7d25ba90c *src/MCMCfcds.h +40962cb2912c4c8de9bc61f5592b38cb *src/MCMCfcds.h db0074b59394f540436a57cbe7d03544 *src/MCMChierBetaBinom.cc 0ee00772334d6d8a2d95b800d64e62d7 *src/MCMChierEI.cc 0863bf6fe4fccee125a23f124f5a3c63 *src/MCMCirtKdHet.cc @@ -156,7 +165,7 @@ 7aaedc47c770e41d64ddde202910cfda *src/MCMCnbutil.cc 48367d601c0f1e1bbe2ef042a1b2024e *src/MCMCnbutil.h fb5a3f2c4ece479ba75dd03d86b3a9e0 *src/MCMCordfactanal.cc -54df73ab50502395c620335181065d97 *src/MCMCpack_init.c +9b63d3a3a4e5e27fcb237e33be8666e8 *src/MCMCpack_init.c ef5566f824887d9bbef5b76f1b64d839 *src/MCMCrng.h f3e2f911f14c895bd0ad382501e6a5e8 *src/Makevars.in 80e927d7e94f5aa4466584d560ca29e6 *src/Makevars.win @@ -181,6 +190,9 @@ 59b5a65221e166974d605639d1bc95a1 *src/cMCMCnegbinChange.cc cb566ffb11e4b158db52466002f021ba *src/cMCMCoprobit.cc 80b1ff187408bd55b7b43febf7dbea25 *src/cMCMCoprobitChange.cc +4470f19fcc67cbb886f67963566721b9 *src/cMCMCpaircompare.cc +425557af9c99c966737cdebeeeee0c59 *src/cMCMCpaircompare2d.cc +ab7bba4e5131d0d72d922512e64c3ae4 *src/cMCMCpaircompare2dDP.cc d3af207d8651d767109691ed48cf0532 *src/cMCMCpoisson.cc e5c2951646f50ada80694d828303adfe *src/cMCMCpoissonChange.cc 23c968d628d3811efb2021b184ce6a32 *src/cMCMCprobit.cc diff -Nru r-cran-mcmcpack-1.5-0/NAMESPACE r-cran-mcmcpack-1.6-0/NAMESPACE --- r-cran-mcmcpack-1.5-0/NAMESPACE 2021-01-19 05:27:01.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/NAMESPACE 2021-10-06 01:57:26.000000000 +0000 @@ -34,6 +34,9 @@ export(MCMCoprobit) export(MCMCoprobitChange) export(MCMCordfactanal) +export(MCMCpaircompare) +export(MCMCpaircompare2d) +export(MCMCpaircompare2dDP) export(MCMCpoisson) export(MCMCpoissonChange) export(MCMCprobit) diff -Nru r-cran-mcmcpack-1.5-0/R/data.R r-cran-mcmcpack-1.6-0/R/data.R --- r-cran-mcmcpack-1.5-0/R/data.R 2019-11-18 19:55:23.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/data.R 2021-07-25 01:53:29.000000000 +0000 @@ -167,3 +167,33 @@ #' #' @keywords datasets NULL + +#' Euro 2016 data +#' +#' Data on head-to-head outcomes from the 2016 UEFA European Football +#' Championship. +#' +#' @name Euro2016 +#' +#' @docType data +#' +#' @format This dataframe contains all of the head-to-head results from +#' Euro 2016. This includes results from both the group stage and the +#' knock-out rounds. +#' \describe{ +#' \item{dummy.rater}{An artificial "dummy" rater equal to 1 for all +#' matches. Included so that \code{Euro2016} can be used directly with +#' \code{MCMCpack}'s models for pairwise comparisons.} +#' \item{team1}{The home team} +#' \item{team2}{The away team } +#' \item{winner}{The winner of the match. \code{NA} if a draw.} +#' } +#' +#' +#' +#' @source \url{https://en.wikipedia.org/wiki/UEFA_Euro_2016} +#' +#' @keywords datasets +NULL + + diff -Nru r-cran-mcmcpack-1.5-0/R/HDPHMMnegbin.R r-cran-mcmcpack-1.6-0/R/HDPHMMnegbin.R --- r-cran-mcmcpack-1.5-0/R/HDPHMMnegbin.R 2021-01-19 05:26:04.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/HDPHMMnegbin.R 2021-10-06 01:44:08.000000000 +0000 @@ -162,7 +162,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/HDPHMMpoisson.R r-cran-mcmcpack-1.6-0/R/HDPHMMpoisson.R --- r-cran-mcmcpack-1.5-0/R/HDPHMMpoisson.R 2021-01-19 05:22:02.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/HDPHMMpoisson.R 2021-10-06 01:44:21.000000000 +0000 @@ -138,7 +138,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/HDPHSMMnegbin.R r-cran-mcmcpack-1.6-0/R/HDPHSMMnegbin.R --- r-cran-mcmcpack-1.5-0/R/HDPHSMMnegbin.R 2021-01-19 05:21:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/HDPHSMMnegbin.R 2021-10-06 01:44:30.000000000 +0000 @@ -180,7 +180,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/hidden.R r-cran-mcmcpack-1.6-0/R/hidden.R --- r-cran-mcmcpack-1.5-0/R/hidden.R 2019-11-18 19:55:23.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/hidden.R 2021-07-25 01:53:29.000000000 +0000 @@ -75,6 +75,61 @@ return( list(Lambda.eq.constraints, Lambda.ineq.constraints, X.names)) } + + + + +## create constraints for pairwise comparison models +"build.pairwise.theta.constraints" <- + function(theta.constraints, cand.codes, n.cand, factors){ + + ## build initial constraint matrices and assign var names + theta.eq.constraints <- matrix(NA, n.cand, factors) + theta.ineq.constraints <- matrix(0, n.cand, factors) + + rownames(theta.eq.constraints) <- cand.codes + rownames(theta.ineq.constraints) <- cand.codes + + ## setup the equality and inequality contraints on theta + if (length(theta.constraints) != 0){ + constraint.names <- names(theta.constraints) + for (i in 1:length(constraint.names)){ + name.i <- constraint.names[i] + theta.constraints.i <- theta.constraints[[i]] + col.index <- theta.constraints.i[[1]] + replace.element <- theta.constraints.i[[2]] + if (is.numeric(replace.element)){ + theta.eq.constraints[rownames(theta.eq.constraints)==name.i, + col.index] <- replace.element + } + if (replace.element=="+"){ + theta.ineq.constraints[rownames(theta.ineq.constraints)==name.i, + col.index] <- 1 + } + if (replace.element=="-"){ + theta.ineq.constraints[rownames(theta.ineq.constraints)==name.i, + col.index] <- -1 + } + } + } + + testmat <- theta.ineq.constraints * theta.eq.constraints + + if (min(is.na(testmat))==0){ + if ( min(testmat[!is.na(testmat)]) < 0){ + cat("Constraints on theta are logically inconsistent.\n") + stop("Please respecify and call ", calling.function(), " again.\n") + } + } + theta.eq.constraints[is.na(theta.eq.constraints)] <- -999 + + return( list(theta.eq.constraints, theta.ineq.constraints)) + } + + + + + # return name of the calling function "calling.function" <- function(parentheses=TRUE) { @@ -716,6 +771,85 @@ return(list(v,S)) } + + + + + +## generate starting values for pairwise comparison thetas +"pairwise.theta.start" <- + function(theta.start, K, factors, theta.eq.constraints, + theta.ineq.constraints){ + + theta <- matrix(0, K, factors) + if (any(is.na(theta.start))){# sets theta to equality constraints & 0s + for (i in 1:K){ + for (j in 1:factors){ + if (theta.eq.constraints[i,j]==-999){ + if(theta.ineq.constraints[i,j]==0){ + theta[i,j] <- 0 + } + if(theta.ineq.constraints[i,j]>0){ + theta[i,j] <- .5 + } + if(theta.ineq.constraints[i,j]<0){ + theta[i,j] <- -.5 + } + } + else theta[i,j] <- theta.eq.constraints[i,j] + } + } + } + else if (is.matrix(theta.start)){ + if (nrow(theta.start)==K && ncol(theta.start)==factors) + theta <- theta.start + else { + cat("theta.start not of correct size for model specification.\n") + stop("Please respecify and call ", calling.function(), " again.\n") + } + for (i in 1:K){ + for (j in 1:factors){ + if (theta.eq.constraints[i,j]==-999){ + if(theta.ineq.constraints[i,j]>0 & theta[i,j] <= 0){ + theta[i,j] <- -1 * theta[i,j] + 0.001 + } + if(theta.ineq.constraints[i,j]<0 & theta[i,j] >= 0){ + theta[i,j] <- -1 * theta[i,j] - 0.001 + } + } + else theta[i,j] <- theta.eq.constraints[i,j] + } + } + } + else if (length(theta.start)==1 && is.numeric(theta.start)){ + theta <- matrix(theta.start, K, factors) + for (i in 1:K){ + for (j in 1:factors){ + if (theta.eq.constraints[i,j]==-999){ + if(theta.ineq.constraints[i,j]>0 & theta[i,j] <= 0){ + theta[i,j] <- -1 * theta[i,j] + 0.001 + } + if(theta.ineq.constraints[i,j]<0 & theta[i,j] >= 0){ + theta[i,j] <- -1 * theta[i,j] - 0.001 + } + } + else theta[i,j] <- theta.eq.constraints[i,j] + } + } + } + else { + cat("theta.start neither NA, matrix, nor scalar.\n") + stop("Please respecify and call ", calling.function, " again.\n") + } + + return(theta) + } + + + + + + # parse formula and return a list that contains the model response # matrix as element one, and the model matrix as element two "parse.formula" <- diff -Nru r-cran-mcmcpack-1.5-0/R/HMMpanelFE.R r-cran-mcmcpack-1.6-0/R/HMMpanelFE.R --- r-cran-mcmcpack-1.5-0/R/HMMpanelFE.R 2021-01-19 02:38:46.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/HMMpanelFE.R 2021-10-06 01:52:53.000000000 +0000 @@ -114,6 +114,10 @@ #' change-point models.'' \emph{Journal of Econometrics}. 86: 221-241. #' #' +#' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. +#' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. +#' #' @keywords models #' #' @examples diff -Nru r-cran-mcmcpack-1.5-0/R/HMMpanelRE.R r-cran-mcmcpack-1.6-0/R/HMMpanelRE.R --- r-cran-mcmcpack-1.5-0/R/HMMpanelRE.R 2021-01-19 04:04:41.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/HMMpanelRE.R 2021-10-06 01:46:21.000000000 +0000 @@ -176,6 +176,10 @@ #' change-point models.'' \emph{Journal of Econometrics}. 86: #' 221-241. #' +#' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. +#' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of +#' Statistical Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. +#' #' @keywords models #' #' @examples diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCbinaryChange.R r-cran-mcmcpack-1.6-0/R/MCMCbinaryChange.R --- r-cran-mcmcpack-1.5-0/R/MCMCbinaryChange.R 2021-01-19 02:41:23.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCbinaryChange.R 2021-10-06 01:50:17.000000000 +0000 @@ -99,7 +99,7 @@ #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of #' Statistical Software}. 42(9): 1-21. -#' \url{https://www.jstatsoft.org/v42/i09/}. +#' \doi{10.18637/jss.v042.i09}. #' #' Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs #' Output.'' \emph{Journal of the American Statistical diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCdynamicEI.R r-cran-mcmcpack-1.6-0/R/MCMCdynamicEI.R --- r-cran-mcmcpack-1.5-0/R/MCMCdynamicEI.R 2021-01-19 02:41:57.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCdynamicEI.R 2021-10-06 01:50:05.000000000 +0000 @@ -133,7 +133,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Radford Neal. 2003. ``Slice Sampling" (with discussion). \emph{Annals of #' Statistics}, 31: 705-767. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCdynamicIRT1d.R r-cran-mcmcpack-1.6-0/R/MCMCdynamicIRT1d.R --- r-cran-mcmcpack-1.5-0/R/MCMCdynamicIRT1d.R 2021-01-19 02:42:48.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCdynamicIRT1d.R 2021-10-06 01:49:53.000000000 +0000 @@ -245,7 +245,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' @keywords models #' diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCfactanal.R r-cran-mcmcpack-1.6-0/R/MCMCfactanal.R --- r-cran-mcmcpack-1.5-0/R/MCMCfactanal.R 2021-01-19 02:43:01.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCfactanal.R 2021-10-06 01:49:41.000000000 +0000 @@ -178,7 +178,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMChierEI.R r-cran-mcmcpack-1.6-0/R/MCMChierEI.R --- r-cran-mcmcpack-1.5-0/R/MCMChierEI.R 2021-01-19 02:43:12.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMChierEI.R 2021-10-06 01:49:33.000000000 +0000 @@ -134,7 +134,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCirt1d.R r-cran-mcmcpack-1.6-0/R/MCMCirt1d.R --- r-cran-mcmcpack-1.5-0/R/MCMCirt1d.R 2021-01-19 02:44:03.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCirt1d.R 2021-10-06 01:49:26.000000000 +0000 @@ -179,7 +179,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCirtHier1d.R r-cran-mcmcpack-1.6-0/R/MCMCirtHier1d.R --- r-cran-mcmcpack-1.5-0/R/MCMCirtHier1d.R 2021-01-19 02:44:15.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCirtHier1d.R 2021-10-06 01:49:15.000000000 +0000 @@ -221,7 +221,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCirtKd.R r-cran-mcmcpack-1.6-0/R/MCMCirtKd.R --- r-cran-mcmcpack-1.5-0/R/MCMCirtKd.R 2021-01-19 02:44:21.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCirtKd.R 2021-10-06 01:49:07.000000000 +0000 @@ -189,7 +189,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCirtKdRob.R r-cran-mcmcpack-1.6-0/R/MCMCirtKdRob.R --- r-cran-mcmcpack-1.5-0/R/MCMCirtKdRob.R 2021-01-19 02:44:30.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCirtKdRob.R 2021-10-06 01:48:59.000000000 +0000 @@ -286,7 +286,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMClogit.R r-cran-mcmcpack-1.6-0/R/MCMClogit.R --- r-cran-mcmcpack-1.5-0/R/MCMClogit.R 2021-01-19 02:44:38.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMClogit.R 2021-10-06 01:48:53.000000000 +0000 @@ -142,7 +142,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCmetrop1R.R r-cran-mcmcpack-1.6-0/R/MCMCmetrop1R.R --- r-cran-mcmcpack-1.5-0/R/MCMCmetrop1R.R 2021-01-19 02:44:43.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCmetrop1R.R 2021-10-06 01:48:45.000000000 +0000 @@ -140,7 +140,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCmixfactanal.R r-cran-mcmcpack-1.6-0/R/MCMCmixfactanal.R --- r-cran-mcmcpack-1.5-0/R/MCMCmixfactanal.R 2021-01-19 02:44:51.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCmixfactanal.R 2021-10-06 01:47:48.000000000 +0000 @@ -220,7 +220,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for #' Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.} diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCmnl.R r-cran-mcmcpack-1.6-0/R/MCMCmnl.R --- r-cran-mcmcpack-1.5-0/R/MCMCmnl.R 2021-01-19 02:45:02.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCmnl.R 2021-10-06 01:47:41.000000000 +0000 @@ -366,7 +366,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCnegbinChange.R r-cran-mcmcpack-1.6-0/R/MCMCnegbinChange.R --- r-cran-mcmcpack-1.5-0/R/MCMCnegbinChange.R 2021-01-19 05:26:31.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCnegbinChange.R 2021-10-06 01:43:53.000000000 +0000 @@ -156,7 +156,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCnegbin.R r-cran-mcmcpack-1.6-0/R/MCMCnegbin.R --- r-cran-mcmcpack-1.5-0/R/MCMCnegbin.R 2021-01-19 02:31:52.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCnegbin.R 2021-10-06 01:52:17.000000000 +0000 @@ -124,7 +124,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCoprobitChange.R r-cran-mcmcpack-1.6-0/R/MCMCoprobitChange.R --- r-cran-mcmcpack-1.5-0/R/MCMCoprobitChange.R 2021-01-19 02:46:04.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCoprobitChange.R 2021-10-06 01:47:25.000000000 +0000 @@ -141,7 +141,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point #' models.'' \emph{Journal of Econometrics}. 86: 221-241. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCoprobit.R r-cran-mcmcpack-1.6-0/R/MCMCoprobit.R --- r-cran-mcmcpack-1.5-0/R/MCMCoprobit.R 2021-01-19 02:45:21.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCoprobit.R 2021-10-06 01:47:32.000000000 +0000 @@ -135,7 +135,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Valen E. Johnson and James H. Albert. 1999. \emph{Ordinal Data Modeling}. #' Springer: New York. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCordfactanal.R r-cran-mcmcpack-1.6-0/R/MCMCordfactanal.R --- r-cran-mcmcpack-1.5-0/R/MCMCordfactanal.R 2021-01-19 02:46:14.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCordfactanal.R 2021-10-06 01:47:18.000000000 +0000 @@ -185,7 +185,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' M. K. Cowles. 1996. ``Accelerating Monte Carlo Markov Chain Convergence for #' Cumulative-link Generalized Linear Models." \emph{Statistics and Computing.} diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCpaircompare2dDP.R r-cran-mcmcpack-1.6-0/R/MCMCpaircompare2dDP.R --- r-cran-mcmcpack-1.5-0/R/MCMCpaircompare2dDP.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCpaircompare2dDP.R 2021-10-06 01:42:11.000000000 +0000 @@ -0,0 +1,751 @@ +########################################################################## +## sample from the posterior distribution of a 2-dimensional pairwise +## comparisons model with Dirichlet process prior in R using +## linked C++ code in Scythe. +## +## Model is: +## +## i = 1,...,n.resp (resondents) +## j = 1,...,n.cand (candidates) +## +## Y_{ijj'} = 1 if i chooses j over j' +## Y_{ijj'} = 0 if i chooses j' over j +## Y_{ijj'} = NA if i chooses neither; +## +## Pr(Y_{ijj'} = 1) = \Phi( z_{i}^{T} [\theta_{j} - \theta_{ j'} ] ) +## z_{i}=[cos(gamma_{i}) sin(gamma_{i})]^{T} +## +## gamma_i has a Dirichlet Process prior. +## We use a finite stick breaking process to approximate the infinite +## Dirichlet Process. +## +## theta_j \overset{ind}{\sim} N_{2}(0, I_{2}) +## (some theta_js truncated above or below 0, or fixed to constants) +## +## +## candidate IDs in columns 2 to 4 need to begin with a letter +## +## This software is distributed under the terms of the GNU GENERAL +## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE +## file for more information. +## +## Based on 1-d code from KQ 3/17/2015 +## Initial 2-d DP code QY 11/14/2019 +## Various improvements and modifications KQ and QY 2019 to 2021 +## Added to MCMCpack KQ 6/26/2021 +## +## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn +## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn, +## and Jong Hee Park +########################################################################## + + +#' Markov Chain Monte Carlo for the Two-Dimensional Pairwise Comparisons +#' Model with Dirichlet Process Prior in Yu and Quinn (2021) +#' +#' This function generates a sample from the posterior distribution of a +#' model for pairwise comparisons data with a probit link. Unlike standard +#' models for pairwise comparisons data, in this model the latent attribute +#' of each item being compared is a vector in two-dimensional Euclidean space. +#' +#' \code{MCMCpaircompare2d} uses the data augmentation approach of Albert and +#' Chib (1993) in conjunction with Gibbs and Metropolis-within-Gibbs steps +#' to fit the model. The user supplies data and a sample from the +#' posterior is returned as an \code{mcmc} object, which can be subsequently +#' analyzed in the \code{coda} package. +#' +#' The simulation is done in compiled C++ code to maximize efficiency. +#' +#' Please consult the \code{coda} package documentation for a comprehensive +#' list of functions that can be used to analyze the posterior sample. +#' +#' The model takes the following form: +#' +#' \deqn{i = 1,...,I \ \ \ \ (raters) } +#' \deqn{j = 1,...,J \ \ \ \ (items) } +#' \deqn{Y_{ijj'} = 1 \ \ if \ \ i \ \ chooses \ \ j \ \ over \ \ j'} +#' \deqn{Y_{ijj'} = 0 \ \ if \ \ i \ \ chooses \ \ j' \ \ over \ \ j} +#' \deqn{Y_{ijj'} = NA \ \ if \ \ i \ \ chooses \ \ neither} +#' +#' \deqn{\Pr(Y_{ijj'} = 1) = \Phi( \mathbf{z}_{i}' [\boldsymbol{\theta}_{j} - +#' \boldsymbol{\theta}_{ j'} ])} +#' \deqn{\mathbf{z}_{i}=[\cos(\gamma_{i}), \ \sin(\gamma_{i})]' } +#' +#' The following priors are assumed: +#' \deqn{\gamma_i \sim G} +#' \deqn{G \sim \mathcal{DP}(\alpha G_0)} +#' \deqn{G_0 = \mathcal{U}nif(0, \pi/2)} +#' \deqn{\alpha \sim \mathcal{G}amma(a_0, b_0)} +#' \deqn{\boldsymbol{\theta}_j \sim +#' \mathcal{N}_{2}(\mathbf{0}, \mathbf{I}_{2})} +#' For identification, some \eqn{\boldsymbol{\theta}_j}s are truncated +#' above or below 0, or fixed to constants. +#' +#' +#' @param pwc.data A data.frame containing the pairwise comparisons data. +#' Each row of \code{pwc.data} corresponds to a single pairwise comparison. +#' \code{pwc.data} needs to have exactly four columns. The first column +#' contains a unique identifier for the rater. Column two contains the unique +#' identifier for the first item being compared. Column three contains the +#' unique identifier for the second item being compared. Column four contains +#' the unique identifier of the item selected from the two items being +#' compared. If a tie occurred, the entry in the fourth column should be NA. +#' \strong{The identifiers in columns 2 through 4 must start with a letter. +#' Examples are provided below.} +#' +#' @param theta.constraints A list specifying possible simple equality or +#' inequality constraints on the item parameters. A +#' typical entry in the list has one of three forms: +#' \code{itemname=list(d,c)} which will constrain the dth dimension of +#' theta for the item named \code{itemname} to be equal to c, +#' \code{itemname=list(d,"+")} which will constrain the dth dimension of +#' theta for the item named \code{itemname} to be positive, and +#' \code{itemname=list(d, "-")} which will constrain the dth dimension of +#' theta for the item named \code{itemname} to be negative. +#' +#' @param burnin The number of burn-in iterations for the sampler. +#' +#' @param mcmc The number of Gibbs iterations for the sampler. +#' +#' @param thin The thinning interval used in the simulation. The number of +#' Gibbs iterations must be divisible by this value. +#' +#' @param verbose A switch which determines whether or not the progress of the +#' sampler is printed to the screen. If \code{verbose} is greater than 0 +#' output is printed to the screen every +#' \code{verbose}th iteration. +#' +#' @param seed The seed for the random number generator. If NA, the Mersenne +#' Twister generator is used with default seed 12345; if an integer is passed +#' it is used to seed the Mersenne twister. The user can also pass a list of +#' length two to use the L'Ecuyer random number generator, which is suitable +#' for parallel computation. The first element of the list is the L'Ecuyer +#' seed, which is a vector of length six or NA (if NA a default seed of +#' \code{rep(12345,6)} is used). The second element of list is a positive +#' substream number. See the MCMCpack specification for more details. +#' +#' @param gamma.start The starting value for the gamma vector. This +#' can either be a scalar or a column vector with dimension equal to the number +#' of raters. If this takes a scalar value, then that value will serve as the +#' starting value for all of the gammas. The default value of NA will set the +#' starting value of each gamma parameter to \eqn{\pi/4}. +#' +#' +#' @param theta.start Starting values for the theta. Can be either a numeric +#' scalar, a J by 2 matrix (where J is the number of items compared), or NA. +#' If a scalar, all theta values are set to that value (except elements already +#' specified via theta.contraints. If NA, then non constrained elements of +#' theta are set equal to 0, elements constrained to be positive are set equal +#' to 0.5, elements constrained to be negative are set equal to -0.5 and +#' elements with equality constraints are set to satisfy those constraints. +#' +#' @param store.theta Should the theta draws be returned? Default is TRUE. +#' +#' @param store.gamma Should the gamma draws be returned? Default is TRUE. +#' +#' @param tune Tuning parameter for the random walk Metropolis proposal for +#' each gamma_i. \code{tune} is the width of the uniform proposal centered at +#' the current value of gamma_i. Must be a positive scalar. +#' +#' @param procrustes Should the theta and gamma draws be post-processed with +#' a Procrustes transformation? Default is FALSE. The Procrustes target matrix +#' is derived from the constrained elements of theta. Each row of theta that +#' has both theta values constrained is part of the of the target matrix. +#' Elements with equality constraints are set to those values. Elements +#' constrained to be positive are set to 1. Elements constrained to be negative +#' are set to -1. If \code{procrustes} is set to \code{TRUE} theta.constraints +#' must be set so that there are at least three rows of theta that have both +#' elements of theta constrained. +#' +#' @param alpha.start The starting value for the DP concentration parameter +#' alpha. Must be a positive scalar. Defaults to 1. If \code{alpha.fixed} is +#' set equal to TRUE, then alpha is held fixed at alpha.start. +#' +#' @param cluster.max The maximum number of clusters allowed in the +#' approximation to the DP prior for gamma. Defaults to 100. Must be a +#' positive integer. +#' +#' @param cluster.mcmc The number of additional MCMC iterations that are done +#' to sample each cluster-specific gamma value within one main MCMC iteration. +#' Must be a positive integer. Defaults to 500. Setting this to a lower value +#' speeds runtime at the cost of (possibly) worse mixing. +#' +#' @param alpha.fixed Logical value indicating whether the DP concentration +#' parameter alpha be held fixed (\code{TRUE}) or estimated (\code{FALSE}). +#' +#' @param a0 The shape parameter of the gamma prior for alpha. This is the +#' same parameterization of the gamma distribution as R's internal +#' \code{rgamma()} function. Only relevant if \code{alpha.fixed} is set equal +#' to \code{FALSE}. Defaults to 1. +#' +#' @param b0 The rate parameter of the gamma prior for alpha. This is the +#' same parameterization of the gamma distribution as R's internal +#' \code{rgamma()} function. Only relevant if \code{alpha.fixed} is set equal +#' to \code{FALSE}. Defaults to 1. +#' +#' @param ... further arguments to be passed +#' +#' +#' @return An mcmc object that contains the posterior sample. This object can +#' be summarized by functions provided by the coda package. Most of the column +#' names of the mcmc object are self explanatory. Note however that the columns +#' with names of the form "cluster.[raterID]" give the cluster membership of +#' each rater at each stored MCMC iteration. Because of the possibility of +#' label switching, the particular values of these cluster membership variables +#' are not meaningful. What is meaningful is whether two raters share the same +#' cluster membership value at a particular MCMC iteration. This indicates +#' that those two raters were clustered together during that iteration. +#' Finally, note that the "n.clusters" column gives the number of distinct +#' gamma values at each iteration, i.e. the number of clusters at that +#' iteration. +#' +#' @author Qiushi Yu and +#' Kevin M. Quinn +#' +#' @seealso \code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, +#' \code{\link[MCMCpack]{MCMCpaircompare}}, +#' \code{\link[MCMCpack]{MCMCpaircompare2dDP}} +#' +#' @references Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary +#' and Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, +#' 669-679 +#' +#' Yu, Qiushi and Kevin M. Quinn. 2021. ``A Multidimensional Pairwise +#' Comparison Model for Heterogeneous Perceptions with an Application to +#' Modeling the Perceived Truthfulness of Public Statements on COVID-19.'' +#' University of Michigan Working Paper. +#' +#' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: +#' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. +#' +#' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe +#' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. +#' +#' Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output +#' Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. +#' \url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. +#' +#' @keywords models +#' +#' @examples +#' \dontrun{ +#' ## a synthetic data example +#' set.seed(123) +#' +#' I <- 65 ## number of raters +#' J <- 50 ## number of items to be compared +#' +#' ## 3 clusters: +#' ## raters 1 to 5 put most weight on dimension 1 +#' ## raters 6 to 10 put most weight on dimension 2 +#' ## raters 11 to I put substantial weight on both dimensions +#' gamma.true <- c(rep(0.05, 5), +#' rep(1.50, 5), +#' rep(0.7, I-10) ) +#' theta1.true <- rnorm(J, m=0, s=1) +#' theta2.true <- rnorm(J, m=0, s=1) +#' theta1.true[1] <- 2 +#' theta2.true[1] <- 2 +#' theta1.true[2] <- -2 +#' theta2.true[2] <- -2 +#' theta1.true[3] <- 2 +#' theta2.true[3] <- -2 +#' +#' +#' +#' n.comparisons <- 125 ## number of pairwise comparisons for each rater +#' +#' ## generate synthetic data according to the assumed model +#' rater.id <- NULL +#' item.1.id <- NULL +#' item.2.id <- NULL +#' choice.id <- NULL +#' for (i in 1:I){ +#' for (c in 1:n.comparisons){ +#' rater.id <- c(rater.id, i+100) +#' item.numbers <- sample(1:J, size=2, replace=FALSE) +#' item.1 <- item.numbers[1] +#' item.2 <- item.numbers[2] +#' item.1.id <- c(item.1.id, item.1) +#' item.2.id <- c(item.2.id, item.2) +#' z <- c(cos(gamma.true[i]), sin(gamma.true[i])) +#' eta <- z[1] * (theta1.true[item.1] - theta1.true[item.2]) + +#' z[2] * (theta2.true[item.1] - theta2.true[item.2]) +#' prob.item.1.chosen <- pnorm(eta) +#' u <- runif(1) +#' if (u <= prob.item.1.chosen){ +#' choice.id <- c(choice.id, item.1) +#' } +#' else{ +#' choice.id <- c(choice.id, item.2) +#' } +#' } +#' } +#' item.1.id <- paste("item", item.1.id+100, sep=".") +#' item.2.id <- paste("item", item.2.id+100, sep=".") +#' choice.id <- paste("item", choice.id+100, sep=".") +#' +#' sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id) +#' +#' +#' ## fit the model (should be run for more than 10500 iterations) +#' posterior <- MCMCpaircompare2dDP(pwc.data=sim.data, +#' theta.constraints=list(item.101=list(1,2), +#' item.101=list(2,2), +#' item.102=list(1,-2), +#' item.102=list(2,-2), +#' item.103=list(1,"+"), +#' item.103=list(2,"-")), +#' verbose=100, +#' burnin=500, mcmc=10000, thin=5, +#' cluster.mcmc=10, +#' store.theta=TRUE, store.gamma=TRUE, +#' tune=0.1) +#' +#' +#' +#' +#' +#' theta1.draws <- posterior[, grep("theta1", colnames(posterior))] +#' theta2.draws <- posterior[, grep("theta2", colnames(posterior))] +#' gamma.draws <- posterior[, grep("gamma", colnames(posterior))] +#' +#' theta1.post.med <- apply(theta1.draws, 2, median) +#' theta2.post.med <- apply(theta2.draws, 2, median) +#' gamma.post.med <- apply(gamma.draws, 2, median) +#' +#' theta1.post.025 <- apply(theta1.draws, 2, quantile, prob=0.025) +#' theta1.post.975 <- apply(theta1.draws, 2, quantile, prob=0.975) +#' theta2.post.025 <- apply(theta2.draws, 2, quantile, prob=0.025) +#' theta2.post.975 <- apply(theta2.draws, 2, quantile, prob=0.975) +#' gamma.post.025 <- apply(gamma.draws, 2, quantile, prob=0.025) +#' gamma.post.975 <- apply(gamma.draws, 2, quantile, prob=0.975) +#' +#' +#' +#' ## compare estimates to truth +#' par(mfrow=c(2,2)) +#' plot(theta1.true, theta1.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=theta1.true, x1=theta1.true, +#' y0=theta1.post.025, y1=theta1.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' plot(theta2.true, theta2.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=theta2.true, x1=theta2.true, +#' y0=theta2.post.025, y1=theta2.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' plot(gamma.true, gamma.post.med, xlim=c(0, 1.6), ylim=c(0, 1.6), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=gamma.true, x1=gamma.true, +#' y0=gamma.post.025, y1=gamma.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' +#' ## plot point estimates +#' plot(theta1.post.med, theta2.post.med, +#' xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' for (i in 1:length(gamma.post.med)){ +#' arrows(x0=0, y0=0, +#' x1=cos(gamma.post.med[i]), +#' y1=sin(gamma.post.med[i]), +#' col=rgb(1,0,0,0.2), len=0.05, lwd=0.5) +#' } +#' +#' } +#' @export + + +"MCMCpaircompare2dDP" <- function(pwc.data, theta.constraints=list(), + burnin=1000, mcmc=20000, thin=1, + verbose=0, seed=NA, + gamma.start=NA, + theta.start=NA, + store.theta=TRUE, + store.gamma=FALSE, + tune=0.3, procrustes=FALSE, + alpha.start=1, cluster.max=100, + cluster.mcmc=500, + alpha.fixed=TRUE, + a0=1, b0=1, + ...){ + + ## checks + check.offset(list(...)) + check.mcmc.parameters(burnin, mcmc, thin) + + if (!(procrustes %in% c(TRUE, FALSE))){ + cat("procrustes cannot take a value other than TRUE or FALSE.\n") + stop("Please check function arguments and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + ## convert all columns to character data + pwc.data[,1] <- as.character(pwc.data[,1]) + pwc.data[,2] <- as.character(pwc.data[,2]) + pwc.data[,3] <- as.character(pwc.data[,3]) + pwc.data[,4] <- as.character(pwc.data[,4]) + + ## check input data + if (ncol(pwc.data) != 4){ + cat("pwc.data must have 4 columns. The specified pwc.data does not have 4 columns.\n") + stop("Please check data and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + for (i in 1:nrow(pwc.data)){ + if (!(pwc.data[i,4] %in% c(NA, pwc.data[i,2], pwc.data[i,3]))){ + cat("pwc.data[", i, ",4] is not in {NA, pwc.data[", i, ",2:3]}.\n", sep="") + stop("Please check data and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + } + + + + ## extract key constants from pwc.data + n <- nrow(pwc.data) + n.resp <- length(unique(pwc.data[,1])) + n.cand <- length(unique( c(pwc.data[,2], pwc.data[,3]))) + + + ## convert pwc.data into purely numeric matrix + resp.codes <- sort(unique(pwc.data[,1])) + cand.codes <- sort(unique( c(pwc.data[,2], pwc.data[,3]) )) + pwc.data.numeric <- matrix(-999, nrow(pwc.data), 4) + for (p in 1:n){ + if (pwc.data[p,1] > 0){ + pwc.data.numeric[p,1] <- which(pwc.data[p,1] == resp.codes) + } + if (pwc.data[p,2] > 0){ + pwc.data.numeric[p,2] <- which(pwc.data[p,2] == cand.codes) + } + if (pwc.data[p,3] > 0){ + pwc.data.numeric[p,3] <- which(pwc.data[p,3] == cand.codes) + } + if (pwc.data[p,4] > 0){ + pwc.data.numeric[p,4] <- which(pwc.data[p,4] == cand.codes) + } + } + + + if (length(cluster.max) > 1){ + cat("cluster.max is not a scalar in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (!is.numeric(cluster.max)){ + cat("cluster.max is non-numeric in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (cluster.max <= 0){ + cat("cluster.max takes a value <= 0 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + if (length(cluster.mcmc) > 1){ + cat("cluster.mcmc is not a scalar in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (!is.numeric(cluster.mcmc)){ + cat("cluster.mcmc is non-numeric in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (cluster.mcmc <= 0){ + cat("cluster.mcmc takes a value <= 0 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + + if (length(a0) > 1){ + cat("a0 is not a scalar in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (!is.numeric(a0)){ + cat("a0 is non-numeric in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (a0 <= 0){ + cat("a0 takes a value <= 0 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + + if (length(b0) > 1){ + cat("b0 is not a scalar in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (!is.numeric(b0)){ + cat("b0 is non-numeric in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (b0 <= 0){ + cat("b0 takes a value <= 0 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + ## set up constraints on theta + holder <- build.pairwise.theta.constraints(theta.constraints, + cand.codes, n.cand, 2) + theta.eq.constraints <- holder[[1]] + theta.ineq.constraints <- holder[[2]] + + ## starting values for theta + theta <- pairwise.theta.start(theta.start, n.cand, 2, + theta.eq.constraints, + theta.ineq.constraints) + + ## starting values for alpha + if (length(alpha.start) > 1){ + cat("alpha.start is not a scalar in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (!is.numeric(alpha.start)){ + cat("alpha.start is non-numeric in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (alpha.start <= 0){ + cat("alpha.start takes a value <= 0 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + + ## starting values for cluster.gamma + cluster.gamma<-gamma.start + + if (all(is.na(gamma.start))){ + cluster.gamma <- gamma.start <- rep(pi/4, cluster.max) + } + if (length(gamma.start) < cluster.max){ + cluster.gamma <- rep(gamma.start, length.out=cluster.max) + } + if (!is.numeric(gamma.start)){ + cat("gamma.start is non-numeric in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (min(gamma.start) < 0){ + cat("gamma.start takes a value < 0 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + if (max(gamma.start) > pi/2){ + cat("gamma.start takes a value > pi/2 in MCMCpaircompare2dDP().\n") + stop("Please check specification and try MCMCpaircompare2dDP() again.\n", + call.=FALSE) + } + + ##starting values for judge membership + judge.cluster.membership<-rep(NA, n.resp) + ##starting values gamma_i, which depends on judge i's cluster membership + gamma<-rep(NA,n.resp) + ## I randomly assign judges to 1:cluster.max clusters with equal likelihood + for(i in 1:n.resp){ + judge.cluster.membership[i]<-sample(1:cluster.max, 1) + gamma[i]<-cluster.gamma[judge.cluster.membership[i]] + } + + + ## define holder for posterior sample + if(store.gamma == FALSE & store.theta == TRUE) { + sample <- matrix(data=0, mcmc/thin, 2*n.cand) + } + else if (store.gamma == TRUE & store.theta == FALSE){ + #we store judge i's gamma, judge i's cluster membership, and the # of unique clusters in one iteration + sample <- matrix(data=0, mcmc/thin, n.resp*2+1) + } + else if (store.gamma == TRUE & store.theta == TRUE){ + sample <- matrix(data=0, mcmc/thin, 2*n.cand + n.resp*2+1) + } + else{ + cat("Error: store.gamma == FALSE & store.theta == FALSE.\n") + stop("Please respecify and call MCMCpaircompare() again.\n", + call.=FALSE) + } + #if alpha.fixed=FALSE, we store alpha posterior sample + if(alpha.fixed==FALSE){ + sample<-cbind(sample,0) + } + + + ## define holder for the acceptance rate of each gamma + gamma_accept_rate<-rep(0, n.resp) + + ## seeds + seeds <- form.seeds(seed) + lecuyer <- seeds[[1]] + seed.array <- seeds[[2]] + lecuyer.stream <- seeds[[3]] + + + ## create theta.sub.target and associated row indicator + theta.eq.const.holder <- theta.eq.constraints + theta.ineq.const.holder <- theta.ineq.constraints + theta.eq.const.holder[theta.eq.const.holder[,1] == -999 & + theta.ineq.const.holder[,1] != 0, 1] <- + theta.ineq.const.holder[theta.eq.const.holder[,1] == -999 & + theta.ineq.const.holder[,1] != 0, 1] + theta.eq.const.holder[theta.eq.const.holder[,2] == -999 & + theta.ineq.const.holder[,2] != 0, 2] <- + theta.ineq.const.holder[theta.eq.const.holder[,2] == -999 & + theta.ineq.const.holder[,2] != 0, 2] + + theta.sub.target.indic <- ((theta.eq.const.holder[,1] != -999) & + (theta.eq.const.holder[,2] != -999)) + + theta.sub.target <- theta.eq.const.holder[theta.sub.target.indic,] + + if (procrustes == TRUE){ + if (nrow(theta.sub.target) < 3){ + cat("Error: procrustes == TRUE but theta.constraints has < 3 rows.\n") + stop("Please respecify and call MCMCpaircompare() again.\n", + call.=FALSE) + } + else{ + theta.eq.constraints <- 0 * theta.eq.constraints - 999 + theta.ineq.constraints <- 0 * theta.ineq.constraints + } + } + + + + ## call C++ code to draw sample + posterior <- .C("cMCMCpaircompare2dDP", + sampledata = as.double(sample), + samplerow = as.integer(nrow(sample)), + samplecol = as.integer(ncol(sample)), + pwc.datanumericdata = as.integer(pwc.data.numeric-1),#minus one because respondents and items are indexed from zero in C++ + pwc.datanumericrow = as.integer(nrow(pwc.data.numeric)), + pwc.datanumericcol = as.integer(ncol(pwc.data.numeric)), + burnin = as.integer(burnin), + mcmc = as.integer(mcmc), + clustermcmc = as.integer(cluster.mcmc), + thin = as.integer(thin), + lecuyer = as.integer(lecuyer), + seedarray = as.integer(seed.array), + lecuyerstream = as.integer(lecuyer.stream), + verbose = as.integer(verbose), + thetastartdata = as.double(theta), + thetastartrow = as.integer(nrow(theta)), + thetastartcol = as.integer(ncol(theta)), + gammastartdata = as.double(gamma), + gammastartrow = as.integer(length(gamma)), + gammastartcol = as.integer(1), + clustergammastartdata=as.double(cluster.gamma), + clustergammastartrow= as.integer(length(cluster.gamma)), + clustergammastartcol= as.integer(1), + judgeclustermembershipstartdata=as.integer(judge.cluster.membership-1), #minus one because clusters are indexed from zero in C++ + judgeclustermembershipstartrow=as.integer(length(judge.cluster.membership)), + judgeclustermembershipstartcol= as.integer(1), + tunevalue=as.double(tune), + thetaeqdata = as.double(theta.eq.constraints), + thetaeqrow = as.integer(nrow(theta.eq.constraints)), + thetaeqcol = as.integer(ncol(theta.eq.constraints)), + thetaineqdata = as.double(theta.ineq.constraints), + thetaineqrow = as.integer(nrow(theta.ineq.constraints)), + thetaineqcol = as.integer(ncol(theta.ineq.constraints)), + storegamma = as.integer(store.gamma), + storetheta = as.integer(store.theta), + gammaacceptrate=as.double(gamma_accept_rate), + alpha=as.double(alpha.start), + clustermax=as.integer(cluster.max), + alphafixed=as.integer(alpha.fixed), + a=as.double(a0), + b=as.double(b0), + PACKAGE="MCMCpack" + ) + + ## undo the C++ indexing by 0 + posterior$pwc.datanumericdata <- posterior$pwc.datanumericdata + 1 + + theta.names <- c(paste("theta1.", cand.codes, sep = ""), paste("theta2.", cand.codes, sep = "")) + ## I store theta's in the following way: + ## if we have 1,2,3,4,5 cand.codes, + ## then theta.names is "theta1.1" "theta1.2" "theta1.3" "theta1.4" "theta1.5" + ## "theta2.1" "theta2.2" "theta2.3" "theta2.4" "theta2.5" + gamma.names <- paste("gamma.", resp.codes, sep = "") + + cluster.names <- paste("cluster.", resp.codes, sep = "") + + ## put together matrix and build MCMC object to return + sample <- matrix(posterior$sampledata, posterior$samplerow, + posterior$samplecol, + byrow=FALSE) + output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin) + + names <- NULL + if(store.theta == TRUE) { + names <- c(names, theta.names) + } + if (store.gamma == TRUE){ + names <- c(names, gamma.names, cluster.names, "n.clusters") + } + if (alpha.fixed == FALSE){ + names <- c(names, "alpha") + } + varnames(output) <- names + + + if (procrustes){ + gammas <- output[, grep("gamma", colnames(output))] + thetas <- output[, grep("theta", colnames(output))] + thetas1 <- thetas[,1:n.cand] + thetas2 <- thetas[,(n.cand+1):ncol(thetas)] + for (iter in 1:nrow(thetas)){ + theta.sub <- cbind(thetas1[iter, theta.sub.target.indic], + thetas2[iter, theta.sub.target.indic]) + procrust.out <- procrustes(theta.sub, theta.sub.target, + translation=FALSE, dilation=FALSE) + R <- procrust.out$R + + theta.mat <- cbind(thetas1[iter,], thetas2[iter,]) + theta.mat <- theta.mat %*% R + thetas1[iter, ] <- theta.mat[,1] + thetas2[iter, ] <- theta.mat[,2] + + z <- cbind(cos(gammas[iter,]), sin(gammas[iter,])) + z <- t(R) %*% t(z) + gammas[iter,] <- atan2(z[2,], z[1,]) + } + # output <- mcmc(data=cbind(thetas1, thetas2, gammas), + # start=burnin+1, end=burnin+mcmc, thin=thin) + output[,1:(2*n.cand)]<-cbind(thetas1, thetas2) + output[,(1+2*n.cand):(n.resp+2*n.cand)]<-gammas + } + + + + attr(output,"title") <- + "MCMCpaircompare2dDP Posterior Sample" + gamma_accept_rate <- posterior$gammaacceptrate + attr(output, "gamma_accept_rate") <- gamma_accept_rate + attr(output, "procrustes") <- procrustes + return(output) + + +} ## end MCMCpaircompare + + + + + diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCpaircompare2d.R r-cran-mcmcpack-1.6-0/R/MCMCpaircompare2d.R --- r-cran-mcmcpack-1.5-0/R/MCMCpaircompare2d.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCpaircompare2d.R 2021-10-06 01:42:34.000000000 +0000 @@ -0,0 +1,602 @@ +########################################################################## +## sample from the posterior distribution of a 2-dimensional pairwise +## comparisons model in R using linked C++ code in Scythe. +## +## Model is: +## +## i = 1,...,n.resp (resondents) +## j = 1,...,n.cand (candidates) +## +## Y_{ijj'} = 1 if i chooses j over j' +## Y_{ijj'} = 0 if i chooses j' over j +## Y_{ijj'} = NA if i chooses neither; +## +## Pr(Y_{ijj'} = 1) = \Phi( z_{i}^{T} [\theta_{j} - \theta_{ j'} ] ) +## z_{i}=[cos(gamma_{i}) sin(gamma_{i})]^{T} +## +## gamma_i \overset{iid}{\sim} Unif(0, \pi/2) +## theta_j \overset{ind}{\sim} N_{2}(0, I_{2}) +## (some theta_js truncated above or below 0, or fixed to constants) +## +## +## candidate IDs in columns 2 to 4 need to begin with a letter +## +## This software is distributed under the terms of the GNU GENERAL +## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE +## file for more information. +## +## Based on 1-d code from KQ 3/17/2015 +## Initial 2-d code QY 10/5/2019 +## Various improvements and modifications KQ and QY 2019 to 2021 +## Added to MCMCpack KQ 6/26/2021 +## +## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn +## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn, +## and Jong Hee Park +########################################################################## + + +#' Markov Chain Monte Carlo for the Two-Dimensional Pairwise Comparisons +#' Model in Yu and Quinn (2021) +#' +#' This function generates a sample from the posterior distribution of a +#' model for pairwise comparisons data with a probit link. Unlike standard +#' models for pairwise comparisons data, in this model the latent attribute +#' of each item being compared is a vector in two-dimensional Euclidean space. +#' +#' \code{MCMCpaircompare2d} uses the data augmentation approach of Albert and +#' Chib (1993) in conjunction with Gibbs and Metropolis-within-Gibbs steps +#' to fit the model. The user supplies data and a sample from the +#' posterior is returned as an \code{mcmc} object, which can be subsequently +#' analyzed in the \code{coda} package. +#' +#' The simulation is done in compiled C++ code to maximize efficiency. +#' +#' Please consult the \code{coda} package documentation for a comprehensive +#' list of functions that can be used to analyze the posterior sample. +#' +#' The model takes the following form: +#' +#' \deqn{i = 1,...,I \ \ \ \ (raters) } +#' \deqn{j = 1,...,J \ \ \ \ (items) } +#' \deqn{Y_{ijj'} = 1 \ \ if \ \ i \ \ chooses \ \ j \ \ over \ \ j'} +#' \deqn{Y_{ijj'} = 0 \ \ if \ \ i \ \ chooses \ \ j' \ \ over \ \ j} +#' \deqn{Y_{ijj'} = NA \ \ if \ \ i \ \ chooses \ \ neither} +#' +#' \deqn{\Pr(Y_{ijj'} = 1) = \Phi( \mathbf{z}_{i}' [\boldsymbol{\theta}_{j} - +#' \boldsymbol{\theta}_{ j'} ])} +#' \deqn{\mathbf{z}_{i}=[\cos(\gamma_{i}), \ \sin(\gamma_{i})]' } +#' +#' The following priors are assumed: +#' \deqn{\gamma_i \sim \mathcal{U}nif(0, \ \pi/2)} +#' \deqn{\boldsymbol{\theta}_j \sim +#' \mathcal{N}_{2}(\mathbf{0}, \mathbf{I}_{2})} +#' For identification, some \eqn{\boldsymbol{\theta}_j}s are truncated +#' above or below 0, or fixed to constants. +#' +#' +#' @param pwc.data A data.frame containing the pairwise comparisons data. +#' Each row of \code{pwc.data} corresponds to a single pairwise comparison. +#' \code{pwc.data} needs to have exactly four columns. The first column +#' contains a unique identifier for the rater. Column two contains the unique +#' identifier for the first item being compared. Column three contains the +#' unique identifier for the second item being compared. Column four contains +#' the unique identifier of the item selected from the two items being +#' compared. If a tie occurred, the entry in the fourth column should be NA. +#' \strong{The identifiers in columns 2 through 4 must start with a letter. +#' Examples are provided below.} +#' +#' @param theta.constraints A list specifying possible simple equality or +#' inequality constraints on the item parameters. A +#' typical entry in the list has one of three forms: +#' \code{itemname=list(d,c)} which will constrain the dth dimension of +#' theta for the item named \code{itemname} to be equal to c, +#' \code{itemname=list(d,"+")} which will constrain the dth dimension of +#' theta for the item named \code{itemname} to be positive, and +#' \code{itemname=list(d, "-")} which will constrain the dth dimension of +#' theta for the item named \code{itemname} to be negative. +#' +#' @param burnin The number of burn-in iterations for the sampler. +#' +#' @param mcmc The number of Gibbs iterations for the sampler. +#' +#' @param thin The thinning interval used in the simulation. The number of +#' Gibbs iterations must be divisible by this value. +#' +#' @param verbose A switch which determines whether or not the progress of the +#' sampler is printed to the screen. If \code{verbose} is greater than 0 +#' output is printed to the screen every +#' \code{verbose}th iteration. +#' +#' @param seed The seed for the random number generator. If NA, the Mersenne +#' Twister generator is used with default seed 12345; if an integer is passed +#' it is used to seed the Mersenne twister. The user can also pass a list of +#' length two to use the L'Ecuyer random number generator, which is suitable +#' for parallel computation. The first element of the list is the L'Ecuyer +#' seed, which is a vector of length six or NA (if NA a default seed of +#' \code{rep(12345,6)} is used). The second element of list is a positive +#' substream number. See the MCMCpack specification for more details. +#' +#' @param gamma.start The starting value for the gamma vector. This +#' can either be a scalar or a column vector with dimension equal to the number +#' of raters. If this takes a scalar value, then that value will serve as the +#' starting value for all of the gammas. The default value of NA will set the +#' starting value of each gamma parameter to \eqn{\pi/4}. +#' +#' +#' @param theta.start Starting values for the theta. Can be either a numeric +#' scalar, a J by 2 matrix (where J is the number of items compared), or NA. +#' If a scalar, all theta values are set to that value (except elements already +#' specified via theta.contraints. If NA, then non constrained elements of +#' theta are set equal to 0, elements constrained to be positive are set equal +#' to 0.5, elements constrained to be negative are set equal to -0.5 and +#' elements with equality constraints are set to satisfy those constraints. +#' +#' @param store.theta Should the theta draws be returned? Default is TRUE. +#' +#' @param store.gamma Should the gamma draws be returned? Default is TRUE. +#' +#' @param tune Tuning parameter for the random walk Metropolis proposal for +#' each gamma_i. \code{tune} is the width of the uniform proposal centered at +#' the current value of gamma_i. Must be a positive scalar. +#' +#' @param procrustes Should the theta and gamma draws be post-processed with +#' a Procrustes transformation? Default is FALSE. The Procrustes target matrix +#' is derived from the constrained elements of theta. Each row of theta that +#' has both theta values constrained is part of the of the target matrix. +#' Elements with equality constraints are set to those values. Elements +#' constrained to be positive are set to 1. Elements constrained to be negative +#' are set to -1. If \code{procrustes} is set to \code{TRUE} theta.constraints +#' must be set so that there are at least three rows of theta that have both +#' elements of theta constrained. +#' +#' @param ... further arguments to be passed +#' +#' +#' @return An mcmc object that contains the posterior sample. This object can +#' be summarized by functions provided by the coda package. +#' +#' @author Qiushi Yu and +#' Kevin M. Quinn +#' +#' @seealso \code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, +#' \code{\link[MCMCpack]{MCMCpaircompare}}, +#' \code{\link[MCMCpack]{MCMCpaircompare2dDP}} +#' +#' @references Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary +#' and Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, +#' 669-679 +#' +#' Yu, Qiushi and Kevin M. Quinn. 2021. ``A Multidimensional Pairwise +#' Comparison Model for Heterogeneous Perceptions with an Application to +#' Modeling the Perceived Truthfulness of Public Statements on COVID-19.'' +#' University of Michigan Working Paper. +#' +#' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: +#' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. +#' +#' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe +#' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. +#' +#' Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output +#' Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. +#' \url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. +#' +#' @keywords models +#' +#' @examples +#' +#' \dontrun{ +#' ## a synthetic data example +#' set.seed(123) +#' +#' I <- 65 ## number of raters +#' J <- 50 ## number of items to be compared +#' +#' +#' ## raters 1 to 5 put most weight on dimension 1 +#' ## raters 6 to 10 put most weight on dimension 2 +#' ## raters 11 to I put substantial weight on both dimensions +#' gamma.true <- c(runif(5, 0, 0.1), +#' runif(5, 1.47, 1.57), +#' runif(I-10, 0.58, 0.98) ) +#' theta1.true <- rnorm(J, m=0, s=1) +#' theta2.true <- rnorm(J, m=0, s=1) +#' theta1.true[1] <- 2 +#' theta2.true[1] <- 2 +#' theta1.true[2] <- -2 +#' theta2.true[2] <- -2 +#' theta1.true[3] <- 2 +#' theta2.true[3] <- -2 +#' +#' +#' +#' n.comparisons <- 125 ## number of pairwise comparisons for each rater +#' +#' ## generate synthetic data according to the assumed model +#' rater.id <- NULL +#' item.1.id <- NULL +#' item.2.id <- NULL +#' choice.id <- NULL +#' for (i in 1:I){ +#' for (c in 1:n.comparisons){ +#' rater.id <- c(rater.id, i+100) +#' item.numbers <- sample(1:J, size=2, replace=FALSE) +#' item.1 <- item.numbers[1] +#' item.2 <- item.numbers[2] +#' item.1.id <- c(item.1.id, item.1) +#' item.2.id <- c(item.2.id, item.2) +#' z <- c(cos(gamma.true[i]), sin(gamma.true[i])) +#' eta <- z[1] * (theta1.true[item.1] - theta1.true[item.2]) + +#' z[2] * (theta2.true[item.1] - theta2.true[item.2]) +#' prob.item.1.chosen <- pnorm(eta) +#' u <- runif(1) +#' if (u <= prob.item.1.chosen){ +#' choice.id <- c(choice.id, item.1) +#' } +#' else{ +#' choice.id <- c(choice.id, item.2) +#' } +#' } +#' } +#' item.1.id <- paste("item", item.1.id+100, sep=".") +#' item.2.id <- paste("item", item.2.id+100, sep=".") +#' choice.id <- paste("item", choice.id+100, sep=".") +#' +#' sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id) +#' +#' +#' ## fit the model +#' posterior <- MCMCpaircompare2d(pwc.data=sim.data, +#' theta.constraints=list(item.101=list(1,2), +#' item.101=list(2,2), +#' item.102=list(1,-2), +#' item.102=list(2,-2), +#' item.103=list(1,"+"), +#' item.103=list(2,"-")), +#' verbose=1000, +#' burnin=500, mcmc=20000, thin=10, +#' store.theta=TRUE, store.gamma=TRUE, tune=0.5) +#' +#' +#' +#' +#' +#' theta1.draws <- posterior[, grep("theta1", colnames(posterior))] +#' theta2.draws <- posterior[, grep("theta2", colnames(posterior))] +#' gamma.draws <- posterior[, grep("gamma", colnames(posterior))] +#' +#' theta1.post.med <- apply(theta1.draws, 2, median) +#' theta2.post.med <- apply(theta2.draws, 2, median) +#' gamma.post.med <- apply(gamma.draws, 2, median) +#' +#' theta1.post.025 <- apply(theta1.draws, 2, quantile, prob=0.025) +#' theta1.post.975 <- apply(theta1.draws, 2, quantile, prob=0.975) +#' theta2.post.025 <- apply(theta2.draws, 2, quantile, prob=0.025) +#' theta2.post.975 <- apply(theta2.draws, 2, quantile, prob=0.975) +#' gamma.post.025 <- apply(gamma.draws, 2, quantile, prob=0.025) +#' gamma.post.975 <- apply(gamma.draws, 2, quantile, prob=0.975) +#' +#' +#' +#' ## compare estimates to truth +#' par(mfrow=c(2,2)) +#' plot(theta1.true, theta1.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=theta1.true, x1=theta1.true, +#' y0=theta1.post.025, y1=theta1.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' plot(theta2.true, theta2.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=theta2.true, x1=theta2.true, +#' y0=theta2.post.025, y1=theta2.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' plot(gamma.true, gamma.post.med, xlim=c(0, 1.6), ylim=c(0, 1.6), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=gamma.true, x1=gamma.true, +#' y0=gamma.post.025, y1=gamma.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' +#' ## plot point estimates +#' plot(theta1.post.med, theta2.post.med, +#' xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' for (i in 1:length(gamma.post.med)){ +#' arrows(x0=0, y0=0, +#' x1=cos(gamma.post.med[i]), +#' y1=sin(gamma.post.med[i]), +#' col=rgb(1,0,0,0.2), len=0.05, lwd=0.5) +#' } +#'} +#' @export + + +"MCMCpaircompare2d" <- function(pwc.data, theta.constraints=list(), + burnin=1000, mcmc=20000, thin=1, + verbose=0, seed=NA, + gamma.start=NA, + theta.start=NA, + store.theta=TRUE, + store.gamma=TRUE, + tune=0.3, procrustes=FALSE, + ...){ + + ## checks + check.offset(list(...)) + check.mcmc.parameters(burnin, mcmc, thin) + + if (!is.numeric(tune)){ + cat("tune must be numeric.\n") + stop("Please check function arguments and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + if (length(tune) > 1){ + cat("tune must be a scalar.\n") + stop("Please check function arguments and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + if (tune <= 0){ + cat("tune must be a positive.\n") + stop("Please check function arguments and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + + + if (!(procrustes %in% c(TRUE, FALSE))){ + cat("procrustes cannot take a value other than TRUE or FALSE.\n") + stop("Please check function arguments and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + + ## convert all columns to character data + pwc.data[,1] <- as.character(pwc.data[,1]) + pwc.data[,2] <- as.character(pwc.data[,2]) + pwc.data[,3] <- as.character(pwc.data[,3]) + pwc.data[,4] <- as.character(pwc.data[,4]) + + ## check input data + if (ncol(pwc.data) != 4){ + cat("pwc.data must have 4 columns. The specified pwc.data does not have 4 columns.\n") + stop("Please check data and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + for (i in 1:nrow(pwc.data)){ + if (!(pwc.data[i,4] %in% c(NA, pwc.data[i,2], pwc.data[i,3]))){ + cat("pwc.data[", i, ",4] is not in {NA, pwc.data[", i, ",2:3]}.\n", sep="") + stop("Please check data and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + } + + + + ## extract key constants from pwc.data + n <- nrow(pwc.data) + n.resp <- length(unique(pwc.data[,1])) + n.cand <- length(unique( c(pwc.data[,2], pwc.data[,3]))) + + + ## convert pwc.data into purely numeric matrix + resp.codes <- sort(unique(pwc.data[,1])) + cand.codes <- sort(unique( c(pwc.data[,2], pwc.data[,3]) )) + pwc.data.numeric <- matrix(-999, nrow(pwc.data), 4) + for (p in 1:n){ + if (pwc.data[p,1] > 0){ + pwc.data.numeric[p,1] <- which(pwc.data[p,1] == resp.codes) + } + if (pwc.data[p,2] > 0){ + pwc.data.numeric[p,2] <- which(pwc.data[p,2] == cand.codes) + } + if (pwc.data[p,3] > 0){ + pwc.data.numeric[p,3] <- which(pwc.data[p,3] == cand.codes) + } + if (pwc.data[p,4] > 0){ + pwc.data.numeric[p,4] <- which(pwc.data[p,4] == cand.codes) + } + } + + + ## set up constraints on theta + holder <- build.pairwise.theta.constraints(theta.constraints, + cand.codes, n.cand, 2) + theta.eq.constraints <- holder[[1]] + theta.ineq.constraints <- holder[[2]] + + ## starting values for theta + theta <- pairwise.theta.start(theta.start, n.cand, 2, + theta.eq.constraints, + theta.ineq.constraints) + + + ## starting values for gamma + gamma <- gamma.start + if (all(is.na(gamma.start))){ + gamma <- gamma.start <- rep(pi/4, n.resp) + } + if (length(gamma.start) < n.resp){ + gamma <- rep(gamma.start, length.out=n.resp) + } + if (!is.numeric(gamma.start)){ + cat("gamma.start is non-numeric in MCMCpaircompare2d().\n") + stop("Please check specification and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + if (min(gamma) < 0){ + cat("gamma.start takes a value < 0 in MCMCpaircompare2d().\n") + stop("Please check specification and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + if (max(gamma) > pi/2){ + cat("gamma.start takes a value > pi/2 in MCMCpaircompare2d().\n") + stop("Please check specification and try MCMCpaircompare2d() again.\n", + call.=FALSE) + } + + + + ## define holder for posterior sample + if(store.gamma == FALSE & store.theta == TRUE) { + sample <- matrix(data=0, mcmc/thin, 2*n.cand) + } + else if (store.gamma == TRUE & store.theta == FALSE){ + sample <- matrix(data=0, mcmc/thin, n.resp) + } + else if (store.gamma == TRUE & store.theta == TRUE){ + sample <- matrix(data=0, mcmc/thin, 2*n.cand + n.resp) + } + else{ + cat("Error: store.gamma == FALSE & store.theta == FALSE.\n") + stop("Please respecify and call MCMCpaircompare() again.\n", + call.=FALSE) + } + + + ## define holder for the acceptance rate of each gamma + gamma_accept_rate<-rep(0,n.resp) + + ## seeds + seeds <- form.seeds(seed) + lecuyer <- seeds[[1]] + seed.array <- seeds[[2]] + lecuyer.stream <- seeds[[3]] + + + ## create theta.sub.target and associated row indicator + theta.eq.const.holder <- theta.eq.constraints + theta.ineq.const.holder <- theta.ineq.constraints + theta.eq.const.holder[theta.eq.const.holder[,1] == -999 & + theta.ineq.const.holder[,1] != 0, 1] <- + theta.ineq.const.holder[theta.eq.const.holder[,1] == -999 & + theta.ineq.const.holder[,1] != 0, 1] + theta.eq.const.holder[theta.eq.const.holder[,2] == -999 & + theta.ineq.const.holder[,2] != 0, 2] <- + theta.ineq.const.holder[theta.eq.const.holder[,2] == -999 & + theta.ineq.const.holder[,2] != 0, 2] + + theta.sub.target.indic <- ((theta.eq.const.holder[,1] != -999) & + (theta.eq.const.holder[,2] != -999)) + + theta.sub.target <- theta.eq.const.holder[theta.sub.target.indic,] + + if (procrustes == TRUE){ + if (nrow(theta.sub.target) < 3){ + cat("Error: procrustes == TRUE but theta.constraints has < 3 rows.\n") + stop("Please respecify and call MCMCpaircompare() again.\n", + call.=FALSE) + } + else{ + theta.eq.constraints <- 0 * theta.eq.constraints - 999 + theta.ineq.constraints <- 0 * theta.ineq.constraints + } + } + + + + ## call C++ code to draw sample + posterior <- .C("cMCMCpaircompare2d", + sampledata = as.double(sample), + samplerow = as.integer(nrow(sample)), + samplecol = as.integer(ncol(sample)), + pwc.datanumericdata = as.integer(pwc.data.numeric-1),#minus one because respondents and items are indexed from zero in C++ + pwc.datanumericrow = as.integer(nrow(pwc.data.numeric)), + pwc.datanumericcol = as.integer(ncol(pwc.data.numeric)), + burnin = as.integer(burnin), + mcmc = as.integer(mcmc), + thin = as.integer(thin), + lecuyer = as.integer(lecuyer), + seedarray = as.integer(seed.array), + lecuyerstream = as.integer(lecuyer.stream), + verbose = as.integer(verbose), + thetastartdata = as.double(theta), + thetastartrow = as.integer(nrow(theta)), + thetastartcol = as.integer(ncol(theta)), + gammastartdata = as.double(gamma), + gammastartrow = as.integer(length(gamma)), + gammastartcol = as.integer(1), + tunevalue=as.double(tune), + thetaeqdata = as.double(theta.eq.constraints), + thetaeqrow = as.integer(nrow(theta.eq.constraints)), + thetaeqcol = as.integer(ncol(theta.eq.constraints)), + thetaineqdata = as.double(theta.ineq.constraints), + thetaineqrow = as.integer(nrow(theta.ineq.constraints)), + thetaineqcol = as.integer(ncol(theta.ineq.constraints)), + storegamma = as.integer(store.gamma), + storetheta = as.integer(store.theta), + gammaacceptrate=as.double(gamma_accept_rate), + PACKAGE="MCMCpack" + ) + + ## undo the C++ indexing by 0 + posterior$pwc.datanumericdata <- posterior$pwc.datanumericdata + 1 + + theta.names <- c(paste("theta1.", cand.codes, sep = ""), paste("theta2.", cand.codes, sep = "")) + ## I store theta's in the following way: + ## if we have 1,2,3,4,5 cand.codes, + ## then theta.names is "theta1.1" "theta1.2" "theta1.3" "theta1.4" "theta1.5" + ## "theta2.1" "theta2.2" "theta2.3" "theta2.4" "theta2.5" + gamma.names <- paste("gamma.", resp.codes, sep = "") + + ## put together matrix and build MCMC object to return + sample <- matrix(posterior$sampledata, posterior$samplerow, + posterior$samplecol, + byrow=FALSE) + output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin) + + names <- NULL + if(store.theta == TRUE) { + names <- c(names, theta.names) + } + if (store.gamma == TRUE){ + names <- c(names, gamma.names) + } + varnames(output) <- names + + + if (procrustes){ + gammas <- output[, grep("gamma", colnames(output))] + thetas <- output[, grep("theta", colnames(output))] + thetas1 <- thetas[,1:n.cand] + thetas2 <- thetas[,(n.cand+1):ncol(thetas)] + for (iter in 1:nrow(thetas)){ + theta.sub <- cbind(thetas1[iter, theta.sub.target.indic], + thetas2[iter, theta.sub.target.indic]) + procrust.out <- procrustes(theta.sub, theta.sub.target, + translation=FALSE, dilation=FALSE) + R <- procrust.out$R + + theta.mat <- cbind(thetas1[iter,], thetas2[iter,]) + theta.mat <- theta.mat %*% R + thetas1[iter, ] <- theta.mat[,1] + thetas2[iter, ] <- theta.mat[,2] + + z <- cbind(cos(gammas[iter,]), sin(gammas[iter,])) + z <- t(R) %*% t(z) + gammas[iter,] <- atan2(z[2,], z[1,]) + } + output <- mcmc(data=cbind(thetas1, thetas2, gammas), + start=burnin+1, end=burnin+mcmc, thin=thin) + } + + + + attr(output,"title") <- + "MCMCpaircompare2d Posterior Sample" + gamma_accept_rate <- posterior$gammaacceptrate + attr(output, "gamma_accept_rate") <- gamma_accept_rate + attr(output, "procrustes") <- procrustes + return(output) + + +} ## end MCMCpaircompare + + + + + diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCpaircompare.R r-cran-mcmcpack-1.6-0/R/MCMCpaircompare.R --- r-cran-mcmcpack-1.5-0/R/MCMCpaircompare.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCpaircompare.R 2021-10-06 01:41:47.000000000 +0000 @@ -0,0 +1,523 @@ +########################################################################## +## sample from the posterior distribution of a pairwise comparisons model +## in R using linked C++ code in Scythe. +## +## Model is: +## +## i = 1,...,n.resp (resondents) +## j = 1,...,n.cand (candidates) +## +## Y_{ijj'} = 1 if i chooses j over j' +## Y_{ijj'} = 0 if i chooses j' over j +## Y_{ijj'} = NA if i chooses neither; +## +## Pr(Y_{ijj'} = 1) = \Phi( \alpha_{i} [\theta_{j} - \theta_{ j'} ] ) +## +## alpha_i \overset{iid}{\sim} N(a, A^{-1}) +## theta_j \overset{ind}{\sim} N(0, 1) +## (some theta_js truncated above or below 0, or fixed to constants) +## +## +## candidate IDs in columns 2 to 4 need to begin with a letter +## +## This software is distributed under the terms of the GNU GENERAL +## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE +## file for more information. +## +## Original code KQ 3/17/2015 +## Added to MCMCpack KQ 6/24/2021 +## +## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn +## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn, +## and Jong Hee Park +########################################################################## + +#' Markov Chain Monte Carlo for a Pairwise Comparisons Model with Probit Link +#' +#' This function generates a sample from the posterior distribution of a +#' model for pairwise comparisons data with a probit link. Thurstone's model +#' is a special case of this model when the \eqn{\alpha} parameter is fixed at +#' 1. +#' +#' \code{MCMCpaircompare} uses the data augmentation approach of Albert and +#' Chib (1993). The user supplies data and priors, and a sample from the +#' posterior is returned as an \code{mcmc} object, which can be subsequently +#' analyzed in the \code{coda} package. +#' +#' The simulation is done in compiled C++ code to maximize efficiency. +#' +#' Please consult the \code{coda} package documentation for a comprehensive +#' list of functions that can be used to analyze the posterior sample. +#' +#' The model takes the following form: +#' +#' \deqn{i = 1,...,I \ \ \ \ (raters) } +#' \deqn{j = 1,...,J \ \ \ \ (items) } +#' \deqn{Y_{ijj'} = 1 \ \ if \ \ i \ \ chooses \ \ j \ \ over \ \ j'} +#' \deqn{Y_{ijj'} = 0 \ \ if \ \ i \ \ chooses \ \ j' \ \ over \ \ j} +#' \deqn{Y_{ijj'} = NA \ \ if \ \ i \ \ chooses \ \ neither} +#' +#' \deqn{Pr(Y_{ijj'} = 1) = \Phi( \alpha_{i} [\theta_{j} - \theta_{ j'} ] ) } +#' +#' The following Gaussian priors are assumed: +#' \deqn{\alpha_i \sim \mathcal{N}(a, A^{-1})} +#' \deqn{\theta_j \sim \mathcal{N}(0, 1)} +#' For identification, some \eqn{\theta_j}s are truncated above or below 0, +#' or fixed to constants. +#' +#' +#' @param pwc.data A data.frame containing the pairwise comparisons data. +#' Each row of \code{pwc.data} corresponds to a single pairwise comparison. +#' \code{pwc.data} needs to have exactly four columns. The first column +#' contains a unique identifier for the rater. Column two contains the unique +#' identifier for the first item being compared. Column three contains the +#' unique identifier for the second item being compared. Column four contains +#' the unique identifier of the item selected from the two items being +#' compared. If a tie occurred, the entry in the fourth column should be NA. +#' For applications without raters (such as sports competitions) all entries +#' in the first column should be set to a single value and \code{alpha.fixed} +#' (see below) should be set to \code{TRUE}. \strong{The identifiers in +#' columns 2 through 4 must start with a letter. Examples are provided below.} +#' +#' @param theta.constraints A list specifying possible simple equality or +#' inequality constraints on the item parameters. A typical entry in the +#' list has one of three forms: \code{itemname=c} which will constrain the +#' item parameter for the item named \code{itemname} to be equal to c, +#' \code{itemname="+"} which will constrain the item parameter for the +#' item named \code{itemname} to be positive, and \code{itemname="-"} which +#' will constrain the item parameter for the item named \code{itemname} to +#' be negative. +#' +#' @param alpha.fixed Should alpha be fixed to a constant value of 1 for all +#' raters? Default is FALSE. If set to FALSE, an alpha value is estimated for +#' each rater. +#' +#' @param burnin The number of burn-in iterations for the sampler. +#' +#' @param mcmc The number of Gibbs iterations for the sampler. +#' +#' @param thin The thinning interval used in the simulation. The number of +#' Gibbs iterations must be divisible by this value. +#' +#' @param verbose A switch which determines whether or not the progress of the +#' sampler is printed to the screen. If \code{verbose} is greater than 0 +#' output is printed to the screen every +#' \code{verbose}th iteration. +#' +#' @param seed The seed for the random number generator. If NA, the Mersenne +#' Twister generator is used with default seed 12345; if an integer is passed +#' it is used to seed the Mersenne twister. The user can also pass a list of +#' length two to use the L'Ecuyer random number generator, which is suitable +#' for parallel computation. The first element of the list is the L'Ecuyer +#' seed, which is a vector of length six or NA (if NA a default seed of +#' \code{rep(12345,6)} is used). The second element of list is a positive +#' substream number. See the MCMCpack specification for more details. +#' +#' @param alpha.start The starting value for the alpha vector. This +#' can either be a scalar or a column vector with dimension equal to the number +#' of alphas. If this takes a scalar value, then that value will serve as the +#' starting value for all of the alphas. The default value of NA will set the +#' starting value of each alpha parameter to 1. +#' +#' @param a The prior mean of alpha. Must be a scalar. Default is 0. +#' +#' @param A The prior precision of alpha. Must be a positive scalar. +#' Default is 0.25 (prior variance is 4). +#' +#' @param store.theta Should the theta draws be returned? Default is TRUE. +#' +#' @param store.alpha Should the alpha draws be returned? Default is FALSE. +#' +#' @param ... further arguments to be passed +#' +#' +#' @return An mcmc object that contains the posterior sample. This object can +#' be summarized by functions provided by the coda package. +#' +#' +#' @seealso \code{\link[coda]{plot.mcmc}},\code{\link[coda]{summary.mcmc}}, +#' \code{\link[MCMCpack]{MCMCpaircompare2d}}, +#' \code{\link[MCMCpack]{MCMCpaircompare2dDP}} +#' +#' @references Albert, J. H. and S. Chib. 1993. ``Bayesian Analysis of Binary +#' and Polychotomous Response Data.'' \emph{J. Amer. Statist. Assoc.} 88, +#' 669-679 +#' +#' Yu, Qiushi and Kevin M. Quinn. 2021. ``A Multidimensional Pairwise +#' Comparison Model for Heterogeneous Perception with an Application to +#' Modeling the Perceived Truthfulness of Public Statements on COVID-19.'' +#' University of Michigan Working Paper. +#' +#' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: +#' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. +#' +#' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe +#' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. +#' +#' Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output +#' Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. +#' \url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. +#' +#' @keywords models +#' +#' @examples +#' +#' \dontrun{ +#' ## Euro 2016 example +#' data(Euro2016) +#' +#' posterior1 <- MCMCpaircompare(pwc.data=Euro2016, +#' theta.constraints=list(Ukraine="-", +#' Portugal="+"), +#' alpha.fixed=TRUE, +#' verbose=10000, +#' burnin=10000, mcmc=500000, thin=100, +#' store.theta=TRUE, store.alpha=FALSE) +#' +#' ## alternative identification constraints +#' posterior2 <- MCMCpaircompare(pwc.data=Euro2016, +#' theta.constraints=list(Ukraine="-", +#' Portugal=1), +#' alpha.fixed=TRUE, +#' verbose=10000, +#' burnin=10000, mcmc=500000, thin=100, +#' store.theta=TRUE, store.alpha=FALSE) +#' +#' +#' +#' +#' +#' +#' +#' +#' ## a synthetic data example with estimated rater-specific parameters +#' set.seed(123) +#' +#' I <- 65 ## number of raters +#' J <- 50 ## number of items to be compared +#' +#' +#' ## raters 1 to 5 have less sensitivity to stimuli than raters 6 through I +#' alpha.true <- c(rnorm(5, m=0.2, s=0.05), rnorm(I - 5, m=1, s=0.1)) +#' theta.true <- sort(rnorm(J, m=0, s=1)) +#' +#' n.comparisons <- 125 ## number of pairwise comparisons for each rater +#' +#' ## generate synthetic data according to the assumed model +#' rater.id <- NULL +#' item.1.id <- NULL +#' item.2.id <- NULL +#' choice.id <- NULL +#' for (i in 1:I){ +#' for (c in 1:n.comparisons){ +#' rater.id <- c(rater.id, i+100) +#' item.numbers <- sample(1:J, size=2, replace=FALSE) +#' item.1 <- item.numbers[1] +#' item.2 <- item.numbers[2] +#' item.1.id <- c(item.1.id, item.1) +#' item.2.id <- c(item.2.id, item.2) +#' eta <- alpha.true[i] * (theta.true[item.1] - theta.true[item.2]) +#' prob.item.1.chosen <- pnorm(eta) +#' u <- runif(1) +#' if (u <= prob.item.1.chosen){ +#' choice.id <- c(choice.id, item.1) +#' } +#' else{ +#' choice.id <- c(choice.id, item.2) +#' } +#' } +#' } +#' item.1.id <- paste("item", item.1.id+100, sep=".") +#' item.2.id <- paste("item", item.2.id+100, sep=".") +#' choice.id <- paste("item", choice.id+100, sep=".") +#' +#' sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id) +#' +#' +#' ## fit the model +#' posterior <- MCMCpaircompare(pwc.data=sim.data, +#' theta.constraints=list(item.101=-2, +#' item.150=2), +#' alpha.fixed=FALSE, +#' verbose=10000, +#' a=0, A=0.5, +#' burnin=10000, mcmc=200000, thin=100, +#' store.theta=TRUE, store.alpha=TRUE) +#' +#' theta.draws <- posterior[, grep("theta", colnames(posterior))] +#' alpha.draws <- posterior[, grep("alpha", colnames(posterior))] +#' +#' theta.post.med <- apply(theta.draws, 2, median) +#' alpha.post.med <- apply(alpha.draws, 2, median) +#' +#' theta.post.025 <- apply(theta.draws, 2, quantile, prob=0.025) +#' theta.post.975 <- apply(theta.draws, 2, quantile, prob=0.975) +#' alpha.post.025 <- apply(alpha.draws, 2, quantile, prob=0.025) +#' alpha.post.975 <- apply(alpha.draws, 2, quantile, prob=0.975) +#' +#' ## compare estimates to truth +#' par(mfrow=c(1,2)) +#' plot(theta.true, theta.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=theta.true, x1=theta.true, +#' y0=theta.post.025, y1=theta.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' plot(alpha.true, alpha.post.med, xlim=c(0, 1.2), ylim=c(0, 3), +#' col=rgb(0,0,0,0.3)) +#' segments(x0=alpha.true, x1=alpha.true, +#' y0=alpha.post.025, y1=alpha.post.975, +#' col=rgb(0,0,0,0.3)) +#' abline(0, 1, col=rgb(1,0,0,0.5)) +#' +#' } +#' +#' @export + + +"MCMCpaircompare" <- function(pwc.data, theta.constraints=list(), + alpha.fixed=FALSE, + burnin=1000, mcmc=20000, thin=1, + verbose=0, seed=NA, + alpha.start=NA, + a=0, A=0.25, + store.theta=TRUE, + store.alpha=FALSE, + ...){ + + ## checks + check.offset(list(...)) + check.mcmc.parameters(burnin, mcmc, thin) + if (!is.logical(alpha.fixed)){ + cat("alpha.fixed must be a logical value.\n") + stop("Please check data and try MCMCpaircompare() again.\n", + call.=FALSE) + } + + ## convert all columns to character data + pwc.data[,1] <- as.character(pwc.data[,1]) + pwc.data[,2] <- as.character(pwc.data[,2]) + pwc.data[,3] <- as.character(pwc.data[,3]) + pwc.data[,4] <- as.character(pwc.data[,4]) + + ## check input data + if (ncol(pwc.data) != 4){ + cat("pwc.data must have 4 columns. The specified pwc.data does not have 4 columns.\n") + stop("Please check data and try MCMCpaircompare() again.\n", + call.=FALSE) + } + for (i in 1:nrow(pwc.data)){ + if (!(pwc.data[i,4] %in% c(NA, pwc.data[i,2], pwc.data[i,3]))){ + cat("pwc.data[", i, ",4] is not in {NA, pwc.data[", i, ",2:3]}.\n", sep="") + stop("Please check data and try MCMCpaircompare() again.\n", + call.=FALSE) + } + } + if (!is.numeric(a) | length(a) !=1){ + cat("a must be a scalar.\n") + stop("Please check specification and try MCMCpaircompare() again.\n", + call.=FALSE) + } + if (!is.numeric(A) | length(A) !=1 | A <= 0){ + cat("A must be a positive scalar.\n") + stop("Please check specification and try MCMCpaircompare() again.\n", + call.=FALSE) + } + + + + ## extract key constants from pwc.data + n <- nrow(pwc.data) + n.resp <- length(unique(pwc.data[,1])) + n.cand <- length(unique( c(pwc.data[,2], pwc.data[,3]))) + + + ## convert pwc.data into purely numeric matrix + resp.codes <- sort(unique(pwc.data[,1])) + cand.codes <- sort(unique( c(pwc.data[,2], pwc.data[,3]) )) + pwc.data.numeric <- matrix(-999, nrow(pwc.data), 4) + for (p in 1:n){ + if (!is.na(pwc.data[p,1])){ + pwc.data.numeric[p,1] <- which(pwc.data[p,1] == resp.codes) + } + if (!is.na(pwc.data[p,2])){ + pwc.data.numeric[p,2] <- which(pwc.data[p,2] == cand.codes) + } + if (!is.na(pwc.data[p,3])){ + pwc.data.numeric[p,3] <- which(pwc.data[p,3] == cand.codes) + } + if (!is.na(pwc.data[p,4])){ + pwc.data.numeric[p,4] <- which(pwc.data[p,4] == cand.codes) + } + } + + + ## set up constraints on theta + if(length(theta.constraints) != 0) { + for (i in 1:length(theta.constraints)){ + theta.constraints[[i]] <- + list(as.integer(1), theta.constraints[[i]][1]) + } + } + theta.eq.constraints <- matrix(NA, n.cand) + theta.ineq.constraints <- matrix(0, n.cand) + rownames(theta.eq.constraints) <- cand.codes + rownames(theta.ineq.constraints) <- cand.codes + if (length(theta.constraints) != 0){ + constraint.names <- names(theta.constraints) + for (i in 1:length(constraint.names)){ + name.i <- constraint.names[i] + theta.constraints.i <- theta.constraints[[i]] + cand.index <- theta.constraints.i[[1]] + replace.element <- theta.constraints.i[[2]] + if (is.numeric(replace.element)){ + theta.eq.constraints[rownames(theta.eq.constraints)==name.i, + cand.index] <- replace.element + } + if (replace.element=="+"){ + theta.ineq.constraints[rownames(theta.ineq.constraints)==name.i, + cand.index] <- 1 + } + if (replace.element=="-"){ + theta.ineq.constraints[rownames(theta.ineq.constraints)==name.i, + cand.index] <- -1 + } + } + } + + testmat <- theta.ineq.constraints * theta.eq.constraints + + if (min(is.na(testmat))==0){ + if ( min(testmat[!is.na(testmat)]) < 0){ + cat("Constraints on theta are logically inconsistent.\n") + stop("Please respecify and call ", calling.function(), " again.\n") + } + } + theta.eq.constraints[is.na(theta.eq.constraints)] <- -999 + + + ## starting values for theta + theta.start <- rep(0, n.cand) + + for (j in 1:n.cand){ + cand.code.j <- cand.codes[j] + if (theta.eq.constraints[cand.code.j,1] != -999){ + theta.start[j] <- theta.eq.constraints[cand.code.j,1] + } + if (theta.ineq.constraints[cand.code.j,1] != 0){ + theta.start[j] <- theta.ineq.constraints[cand.code.j,1] + } + } + + + ## starting values for alpha + if (is.na(alpha.start)){ + alpha.start <- rep(1, n.resp) + } + if (length(alpha.start) < n.resp){ + alpha.start <- rep(alpha.start, length.out=n.resp) + } + if (!is.numeric(alpha.start)){ + cat("alpha.start is non-numeric in MCMCpaircompare().\n") + stop("Please check specification and try MCMCpaircompare() again.\n", + call.=FALSE) + } + + if (alpha.fixed){ + alpha.start <- rep(1, n.resp) + } + + + ## define holder for posterior sample + if(store.alpha == FALSE & store.theta == TRUE) { + sample <- matrix(data=0, mcmc/thin, n.cand) + } + else if (store.alpha == TRUE & store.theta == FALSE){ + sample <- matrix(data=0, mcmc/thin, n.resp) + } + else if (store.alpha == TRUE & store.theta == TRUE){ + sample <- matrix(data=0, mcmc/thin, n.cand + n.resp) + } + else{ + cat("Error: store.alpha == FALSE & store.theta == FALSE.\n") + stop("Please respecify and call MCMCpaircompare() again.\n", + call.=FALSE) + } + + + ## seeds + seeds <- form.seeds(seed) + lecuyer <- seeds[[1]] + seed.array <- seeds[[2]] + lecuyer.stream <- seeds[[3]] + + + + ## call C++ code to draw sample + posterior <- .C("cMCMCpaircompare", + sampledata = as.double(sample), + samplerow = as.integer(nrow(sample)), + samplecol = as.integer(ncol(sample)), + pwc.datanumericdata = as.integer(pwc.data.numeric-1), + pwc.datanumericrow = as.integer(nrow(pwc.data.numeric)), + pwc.datanumericcol = as.integer(ncol(pwc.data.numeric)), + alphafixed = as.integer(alpha.fixed), + burnin = as.integer(burnin), + mcmc = as.integer(mcmc), + thin = as.integer(thin), + lecuyer = as.integer(lecuyer), + seedarray = as.integer(seed.array), + lecuyerstream = as.integer(lecuyer.stream), + verbose = as.integer(verbose), + thetastartdata = as.double(theta.start), + thetastartrow = as.integer(length(theta.start)), + thetastartcol = as.integer(1), + astartdata = as.double(alpha.start), + astartrow = as.integer(length(alpha.start)), + astartcol = as.integer(1), + a=as.double(a), + A=as.double(A), + thetaeqdata = as.double(theta.eq.constraints), + thetaeqrow = as.integer(nrow(theta.eq.constraints)), + thetaeqcol = as.integer(ncol(theta.eq.constraints)), + thetaineqdata = as.double(theta.ineq.constraints), + thetaineqrow = as.integer(nrow(theta.ineq.constraints)), + thetaineqcol = as.integer(ncol(theta.ineq.constraints)), + storealpha = as.integer(store.alpha), + storetheta = as.integer(store.theta), + PACKAGE="MCMCpack" + ) + + ## undo the C++ indexing by 0 + posterior$pwc.datanumericdata <- posterior$pwc.datanumericdata + 1 + + theta.names <- paste("theta.", cand.codes, sep = "") + alpha.names <- paste("alpha.", resp.codes, sep = "") + + ## put together matrix and build MCMC object to return + sample <- matrix(posterior$sampledata, posterior$samplerow, + posterior$samplecol, + byrow=FALSE) + output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin) + + names <- NULL + if(store.theta == TRUE) { + names <- c(names, theta.names) + } + if (store.alpha == TRUE){ + names <- c(names, alpha.names) + } + varnames(output) <- names + attr(output,"title") <- + "MCMCpaircompare Posterior Sample" + return(output) + + +} ## end MCMCpaircompare + + + + + diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCpoissonChange.R r-cran-mcmcpack-1.6-0/R/MCMCpoissonChange.R --- r-cran-mcmcpack-1.5-0/R/MCMCpoissonChange.R 2021-01-19 02:49:08.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCpoissonChange.R 2021-10-06 01:47:01.000000000 +0000 @@ -138,7 +138,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' #' \emph{Journal of the American Statistical Association}. 90: 1313-1321. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCpoisson.R r-cran-mcmcpack-1.6-0/R/MCMCpoisson.R --- r-cran-mcmcpack-1.5-0/R/MCMCpoisson.R 2021-01-19 02:46:24.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCpoisson.R 2021-10-06 01:47:12.000000000 +0000 @@ -116,7 +116,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCprobitChange.R r-cran-mcmcpack-1.6-0/R/MCMCprobitChange.R --- r-cran-mcmcpack-1.5-0/R/MCMCprobitChange.R 2021-01-19 04:03:08.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCprobitChange.R 2021-10-06 01:46:46.000000000 +0000 @@ -120,7 +120,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point #' models.'' \emph{Journal of Econometrics}. 86: 221-241. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCprobit.R r-cran-mcmcpack-1.6-0/R/MCMCprobit.R --- r-cran-mcmcpack-1.5-0/R/MCMCprobit.R 2021-01-19 02:49:30.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCprobit.R 2021-10-06 01:46:53.000000000 +0000 @@ -119,7 +119,7 @@ #' #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. ``MCMCpack: #' Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}. -#' 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' #' \emph{Journal of the American Statistical Association}. 90: 1313-1321. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCregressChange.R r-cran-mcmcpack-1.6-0/R/MCMCregressChange.R --- r-cran-mcmcpack-1.5-0/R/MCMCregressChange.R 2021-01-19 04:04:53.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCregressChange.R 2021-10-06 01:45:47.000000000 +0000 @@ -165,6 +165,10 @@ #' models." \emph{Journal of Econometrics}. 86: 221-241. #' #' +#' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. +#' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of +#' Statistical Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. +#' #' @keywords models #' #' @examples diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCregress.R r-cran-mcmcpack-1.6-0/R/MCMCregress.R --- r-cran-mcmcpack-1.5-0/R/MCMCregress.R 2021-01-19 04:03:26.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCregress.R 2021-10-06 01:46:38.000000000 +0000 @@ -136,7 +136,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Siddhartha Chib. 1995. ``Marginal Likelihood from the Gibbs Output.'' #' \emph{Journal of the American Statistical Association}. 90: 1313-1321. diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCSVDreg.R r-cran-mcmcpack-1.6-0/R/MCMCSVDreg.R --- r-cran-mcmcpack-1.5-0/R/MCMCSVDreg.R 2021-01-19 04:05:37.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCSVDreg.R 2021-10-06 01:45:03.000000000 +0000 @@ -180,7 +180,7 @@ #' Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of #' Statistical Software}. 42(9): 1-21. -#' \url{https://www.jstatsoft.org/v42/i09/}. +#' \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. #' \emph{Scythe Statistical Library 1.0.} diff -Nru r-cran-mcmcpack-1.5-0/R/MCMCtobit.R r-cran-mcmcpack-1.6-0/R/MCMCtobit.R --- r-cran-mcmcpack-1.5-0/R/MCMCtobit.R 2021-01-19 10:43:29.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/R/MCMCtobit.R 2021-10-06 01:43:36.000000000 +0000 @@ -136,7 +136,7 @@ #' #' @references Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011. #' ``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical -#' Software}. 42(9): 1-21. \url{https://www.jstatsoft.org/v42/i09/}. +#' Software}. 42(9): 1-21. \doi{10.18637/jss.v042.i09}. #' #' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe #' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. diff -Nru r-cran-mcmcpack-1.5-0/src/cMCMCpaircompare2d.cc r-cran-mcmcpack-1.6-0/src/cMCMCpaircompare2d.cc --- r-cran-mcmcpack-1.5-0/src/cMCMCpaircompare2d.cc 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/src/cMCMCpaircompare2d.cc 2021-07-25 01:53:29.000000000 +0000 @@ -0,0 +1,569 @@ +////////////////////////////////////////////////////////////////////////// +// MCMCpaircompare.cc is C++ code to estimate a pairwise comparison model. +// +// KQ 3/18/2015 +// +// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn +// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn, +// and Jong Hee Park +////////////////////////////////////////////////////////////////////////// + +#ifndef MCMCPAIRCOMPARE2D_CC +#define MCMCPAIRCOMPARE2D_CC + +#include + +#include "MCMCrng.h" +#include "MCMCfcds.h" +#include "matrix.h" +#include "distributions.h" +#include "stat.h" +#include "la.h" +#include "ide.h" +#include "smath.h" +#include "rng.h" +#include "mersenne.h" +#include "lecuyer.h" + +#include // needed to use Rprintf() +#include // needed to allow user interrupts + +using namespace std; +using namespace scythe; + +/* MCMCpaircompare2d implementation. */ + +// update latent Ystar values in 2-d paired comparisons model +template +void paircompare2d_Ystar_update (Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& theta, const Matrix<>& gamma, + rng& stream){ + + const unsigned int N = MD.rows(); + for (unsigned int i = 0; i < N; ++i){ + const double gamma_i = gamma(MD(i,0)); + const double mean = std::cos(gamma_i)*theta(MD(i, 1), 0) + std::sin(gamma_i)*theta(MD(i, 1), 1) - + std::cos(gamma_i)*theta(MD(i, 2), 0) - std::sin(gamma_i)*theta(MD(i, 2), 1); + + const bool above = MD(i,1) == MD(i,3); // cand 1 chosen + const bool below = MD(i,2) == MD(i,3); // cand 2 chosen + if (above){ + Ystar(i) = stream.rtbnorm_combo(mean, 1.0, 0.0); + } + else if (below){ + Ystar(i) = stream.rtanorm_combo(mean, 1.0, 0.0); + } + else{ + Ystar(i) = stream.rnorm(mean, 1.0); + } + } +} + + + + + +// update theta values in 2-d paired comparisons model +template +void paircompare2d_theta_update (Matrix<>& theta, const Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& gamma, + const Matrix& theta_n, + const Matrix<>& theta_eq, + const Matrix<>& theta_ineq, + const vector< vector < double* > >& theta_Ystar_ptr, + const vector< vector < double* > >& theta_gamma_ptr, + const vector< vector < vector < double* > > >& theta_theta_ptr, + const vector< vector < double > >& theta_sign, + rng& stream){ + + const unsigned int J = theta.rows(); + const Matrix<> I = eye(2); + for (unsigned int j = 0; j < J; ++j){ + Matrix<> X(theta_n[j],2); + Matrix<> z(theta_n[j],1); + + // case 1: no equality constraints at all + if (theta_eq(j,0) == -999 && theta_eq(j,1) == -999){ + for (unsigned int i = 0; i < theta_n[j]; ++i){ + X(i,0)= std::cos(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + X(i,1)= std::sin(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + z(i,0) = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + (std::cos(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][0]+ + std::sin(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][1]); + } + Matrix<> v = invpd(crossprod(X) + I); + Matrix<> m = v * t(X)*z; + Matrix<> v_C = cholesky(v); + // case 1a: no inequality constraints at all + if (theta_ineq(j,0) == 0 && theta_ineq(j,1) == 0){ + Matrix<> theta_j_free = gaxpy(v_C, stream.rnorm(2, 1, 0, 1),m); + theta(j, _) = t(theta_j_free); + } + else{ // case 1b: some inequality constraint holds + Matrix<> theta_j_free = gaxpy(v_C, stream.rnorm(2, 1, 0, 1),m); + double ineq_holds = 0.0; + for (unsigned int jj = 0; jj < 2; ++jj){ + ineq_holds = std::min(ineq_holds, + theta_ineq(j,jj) * theta_j_free(jj)); + } + while (ineq_holds < 0){ + theta_j_free = gaxpy(v_C, stream.rnorm(2, 1, 0, 1), m); + double test = 0; + for (unsigned int jj = 0; jj < 2; ++jj) { + test = std::min(test, theta_ineq(j,jj) * theta_j_free(jj)); + } + ineq_holds = test; + } + theta(j,_) = t(theta_j_free); + } + } + else if (theta_eq(j,0) != -999 && theta_eq(j,1) == -999){ + // case 2: equality constraint on theta_j1 but not on theta_j2 + double xx = 0.0; + double xz = 0.0; + for (unsigned int i = 0; i < theta_n[j]; ++i){ + const double x_ji = std::sin(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + const double z_ji = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + (std::cos(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][0]+ + std::sin(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][1]) + -1*theta_sign[j][i] * theta_eq(j,0) * std::cos(*theta_gamma_ptr[j][i]); + xx += x_ji * x_ji; + xz += x_ji * z_ji; + } + const double v = 1.0/(xx + 1.0); + const double m = v * (xz); + // case 2a: no inequality constraint on theta_j2 + if (theta_ineq(j,1) == 0){ + theta(j,1) = stream.rnorm(m, std::sqrt(v)); + } + else{ // case 2b: inequality constraint on theta_j2 + if (theta_ineq(j,1) > 0) { // theta_j2 > 0 + theta(j,1) = stream.rtbnorm_combo(m, v, 0); + } else { // theta_j2 < 0 + theta(j,1) = stream.rtanorm_combo(m, v, 0); + } + } + } + else if (theta_eq(j,0) == -999 && theta_eq(j,1) != -999){ + // case 3: equality constraint on theta_j2 but not on theta_j1 + double xx = 0.0; + double xz = 0.0; + for (unsigned int i = 0; i < theta_n[j]; ++i){ + const double x_ji = std::cos(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + const double z_ji = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + (std::cos(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][0]+ + std::sin(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][1]) + -1*theta_sign[j][i] * theta_eq(j,1) * std::sin(*theta_gamma_ptr[j][i]); + xx += x_ji * x_ji; + xz += x_ji * z_ji; + } + const double v = 1.0/(xx + 1.0); + const double m = v * (xz); + // case 3a: no inequality constraint on theta_j1 + if (theta_ineq(j,0) == 0){ + theta(j,0) = stream.rnorm(m, std::sqrt(v)); + } + else{ // case 3b: inequality constraint on theta_j1 + if (theta_ineq(j,0) > 0) { // theta_j1 > 0 + theta(j,0) = stream.rtbnorm_combo(m, v, 0); + } else { // theta_j1 < 0 + theta(j,0) = stream.rtanorm_combo(m, v, 0); + } + } + } + else if (theta_eq(j,0) != -999 && theta_eq(j,1) != -999){ + // case 4: equality constraints on both theta_j1 and theta_j2 + theta(j,0) = theta_eq(j,0); + theta(j,1) = theta_eq(j,1); + } + + } // end j loop + +} + + +// update gamma values in 2-d paired comparisons model +template +void paircompare2d_gamma_update (Matrix<>& gamma, + const Matrix& gamma_n, + const vector< vector < double* > >& gamma_Ystar_ptr, + const vector< vector < vector< double* > > >& gamma_theta1_ptr, + const vector< vector < vector< double* > > >& gamma_theta2_ptr, + // const vector< vector < int > >& gamma_choice, + const double& tune, // A random walk step in the MH sampler : step \sim runif(-tune, tune). + vector& gamma_trial, + vector& gamma_accept, + rng& stream){ + + const unsigned int I = gamma.rows(); + for (unsigned int i = 0; i < I; ++i){ + double gamma_i=gamma(i); + double gamma_i_new=gamma(i)+(-2*stream.runif()+1)*tune; + //I use 1.5707959999999 to denote PI/2 + while(gamma_i_new<0 || gamma_i_new>1.5707959999999){ + gamma_i_new=gamma(i)+(-2*stream.runif()+1)*tune; + } + //I assume a uniform prior for gamma, so I only consider likelihood for MH rejection rate. + double log_lik_old=0; + double log_lik_new=0; + double eta_i_j=0;//eta_i_j is the linear function value within the Probit link for judge i's j'th comparison. + double eta_i_j_new=0; + + for (unsigned int j = 0; j < gamma_n[i]; ++j){ + eta_i_j=std::cos(gamma_i)* *gamma_theta1_ptr[i][j][0] + std::sin(gamma_i)**gamma_theta1_ptr[i][j][1] - + std::cos(gamma_i)* *gamma_theta2_ptr[i][j][0] - std::sin(gamma_i)* *gamma_theta2_ptr[i][j][1] ; + log_lik_old += lndnorm(*gamma_Ystar_ptr[i][j], eta_i_j, 1);//lndnorm is the log scale normal density + + eta_i_j_new=std::cos(gamma_i_new)* *gamma_theta1_ptr[i][j][0] + std::sin(gamma_i_new)**gamma_theta1_ptr[i][j][1] - + std::cos(gamma_i_new)* *gamma_theta2_ptr[i][j][0] - std::sin(gamma_i_new)* *gamma_theta2_ptr[i][j][1] ; + log_lik_new += lndnorm(*gamma_Ystar_ptr[i][j], eta_i_j_new, 1); + } + gamma_trial[i]+=1;//A new value has been proposed for gamma i + if(stream.runif() +void MCMCpaircompare2d_impl (rng& stream, + const Matrix& MD, + Matrix<>& theta, Matrix<>& gamma, + const Matrix<>& theta_eq, + const Matrix<>& theta_ineq, + const double tune, + const unsigned int burnin, + const unsigned int mcmc, + const unsigned int thin, + const unsigned int verbose, + const bool storegamma, const bool storetheta, + double* sampledata, + const unsigned int samplesize, + double* gammaacceptrate){ + + + // constants + const unsigned int N = MD.rows();//data matrix + const unsigned int J = theta.rows(); + const unsigned int I = gamma.rows(); + const unsigned int tot_iter = burnin + mcmc; + const unsigned int nsamp = mcmc / thin; + + + // starting values for Ystar + Matrix<> Ystar(N,1); + + + // pre-compute what can be pre-computed + Matrix theta_n(J,1,true,0); // vector of 0s. Item j's total instances of being compared. + Matrix gamma_n(I,1,true,0); // vector of 0s. Judge i's total instances of making comparisons. + for (unsigned int i = 0; i < N; ++i){ + gamma_n(MD(i,0)) += 1; + theta_n(MD(i,1)) += 1; + theta_n(MD(i,2)) += 1; + } + + vector< vector< double* > > theta_Ystar_ptr; + vector< vector< double* > > theta_gamma_ptr; + vector< vector< vector< double* > > > theta_theta_ptr;//the inner vector has two pointers for theta1 and theta2 + vector< vector< double > > theta_sign; + + vector< vector< double* > > gamma_Ystar_ptr; + vector< vector< vector< double* > > > gamma_theta1_ptr; //the inner vector has two pointers for theta1 and theta2 + vector< vector< vector< double* > > > gamma_theta2_ptr; //the inner vector has two pointers for theta1 and theta2 + vector< vector< int > > gamma_choice; + // One inner vector of the above vector of vectors stores judge i's choice between all the pairs he works on. + // 1 indicates picking c1, 0 indicates picking c2 + vector gamma_trial; + vector gamma_accept; + + theta_Ystar_ptr.reserve(J); + theta_gamma_ptr.reserve(J); + theta_theta_ptr.reserve(J); + theta_sign.reserve(J); + + gamma_Ystar_ptr.reserve(I); + gamma_theta1_ptr.reserve(I); + gamma_theta2_ptr.reserve(I); + // gamma_choice.reserve(I); + gamma_trial.reserve(I); + gamma_accept.reserve(I); + for(unsigned int i=0; i Ystar_j_ptr; + vector< double* > gamma_j_ptr; + vector< vector > theta_j_ptr; + vector< double > sign_j; + + Ystar_j_ptr.reserve(theta_n(j)); + gamma_j_ptr.reserve(theta_n(j)); + theta_j_ptr.reserve(theta_n(j)); + sign_j.reserve(theta_n(j)); + /* + Allocate proper length for the empty vectors (Ystar_j_ptr, gamma_j_ptr, theta_j_ptr, sign_j) + Then we save these pre-allocated vectors inside the outer layer vectors. + */ + theta_Ystar_ptr.push_back(Ystar_j_ptr); + theta_gamma_ptr.push_back(gamma_j_ptr); + theta_theta_ptr.push_back(theta_j_ptr); + theta_sign.push_back(sign_j); + } + + for (unsigned int i = 0; i < I; ++i){ + vector< double* > Ystar_i_ptr; + vector< vector< double* > > theta1_i_ptr; + vector< vector< double* > > theta2_i_ptr; + // vector< int > choice_i; + + Ystar_i_ptr.reserve(gamma_n(i)); + theta1_i_ptr.reserve(gamma_n(i)); + theta2_i_ptr.reserve(gamma_n(i)); + // choice_i.reserve(gamma_n(i)); + /* + Allocate proper length for the empty vectors (Ystar_j_ptr, theta1_i_ptr, theta2_i_ptr) + Then we save these pre-allocated vectors inside the outer layer vectors. + */ + gamma_Ystar_ptr.push_back(Ystar_i_ptr); + gamma_theta1_ptr.push_back(theta1_i_ptr); + gamma_theta2_ptr.push_back(theta2_i_ptr); + // gamma_choice.push_back(choice_i); + } + + + for (unsigned int i = 0; i < N; ++i){ + unsigned int resp = MD(i,0); + unsigned int c1 = MD(i,1); + unsigned int c2 = MD(i,2); + /* + Assign the starting values of Ystar, gamma, theta and sign to the above vectors. + This is convenient for accessing wanted objects in later updation. + */ + theta_Ystar_ptr[c1].push_back(&Ystar(i,0)); + theta_gamma_ptr[c1].push_back(&gamma(MD(i,0))); + vector theta_temp1; + theta_temp1.push_back(&theta(MD(i,2),0)); + theta_temp1.push_back(&theta(MD(i,2),1)); + theta_theta_ptr[c1].push_back(theta_temp1); + theta_sign[c1].push_back(1.0); + + theta_Ystar_ptr[c2].push_back(&Ystar(i,0)); + theta_gamma_ptr[c2].push_back(&gamma(MD(i,0))); + vector theta_temp2; + theta_temp2.push_back(&theta(MD(i,1),0)); + theta_temp2.push_back(&theta(MD(i,1),1)); + theta_theta_ptr[c2].push_back(theta_temp2); + theta_sign[c2].push_back(-1.0); + + gamma_Ystar_ptr[resp].push_back(&Ystar(i,0)); + gamma_theta1_ptr[resp].push_back(theta_temp2);//theta_temp2 represents c1 + gamma_theta2_ptr[resp].push_back(theta_temp1);//theta_temp1 represents c2 + + + // const bool above_i = MD(i,1) == MD(i,3); // cand 1 chosen + // const bool below_i = MD(i,2) == MD(i,3); // cand 2 chosen + // if (above_i){ + // gamma_choice[resp].push_back(1); + // } + // else if (below_i){ + // gamma_choice[resp].push_back(0); + // } + + + theta_temp1.clear(); + theta_temp2.clear(); + } + + + + + + // storage matrices (col major order) + Matrix<> theta_store; + Matrix<> gamma_store; + if (storetheta) + theta_store = Matrix<>(nsamp, J*2); + + if (storegamma) + gamma_store = Matrix<>(nsamp, I); + + + + unsigned int count = 0; + // MCMC sampling occurs in this for loop + for (unsigned int iter = 0; iter < tot_iter; ++iter){ + + // The following three functions are defined in MCMCfcds.h + // sample Ystar + paircompare2d_Ystar_update(Ystar, MD, theta, gamma, stream); + + // sample gamma + paircompare2d_gamma_update(gamma, + gamma_n, + gamma_Ystar_ptr, + gamma_theta1_ptr, + gamma_theta2_ptr, + //gamma_choice, + tune, + gamma_trial, + gamma_accept, + stream); + + + // sample theta + paircompare2d_theta_update(theta, Ystar, MD, gamma, theta_n, theta_eq, + theta_ineq, + theta_Ystar_ptr, + theta_gamma_ptr, + theta_theta_ptr, + theta_sign, + stream); + + + // print results to screen + if (verbose > 0 && iter % verbose == 0) { + Rprintf("\n\nMCMCpaircompare2d iteration %i of %i \n", (iter+1), tot_iter); + Rprintf("\nSummary of acceptance rates for gamma:\n"); + double gamma_accept_min = 1.0; + double gamma_accept_max = 0.0; + double gamma_accept_num = 0.0; + double gamma_accept_denom = 0.0; + for (unsigned int i=0; i gamma_accept_max) ? local_accrate : gamma_accept_max; + gamma_accept_num += local_accrate; + gamma_accept_denom += 1.0; + } + Rprintf("Minimum gamma acceptance rate: %7.3f", gamma_accept_min, " "); + Rprintf("\nMean gamma acceptance rate: %7.3f", gamma_accept_num/gamma_accept_denom, " "); + Rprintf("\nMaximum gamma acceptance rate: %7.3f", gamma_accept_max, " "); + Rprintf("\n\n"); + + } + + // store results + if ((iter >= burnin) && ((iter % thin == 0))) { + + // store theta + if (storetheta) { + for (unsigned int j = 0; j < J; ++j) { + theta_store(count, j) = theta(j,0); + theta_store(count, J+j) = theta(j, 1); + } + } + + + // store gamma + if (storegamma) + gamma_store(count, _) = gamma; + + count++; + } + + R_CheckUserInterrupt(); // allow user interrupts + + } // end MCMC sampling loop + + + + + // put output back into sampledata + Matrix<> output; + if(storetheta && ! storegamma) { + output = theta_store; + } else if (storegamma && ! storetheta){ + output = gamma_store; + } else { + output = cbind(theta_store, gamma_store); + } + + for (unsigned int i = 0; i < samplesize; ++i) + sampledata[i] = output[i]; + + for (unsigned int i = 0; i < I; ++i) + gammaacceptrate[i] = gamma_accept[i]/gamma_trial[i]; + +} // end MCMCpaircompare implementation + + + + + + + + + +extern "C" { + + void + cMCMCpaircompare2d(double* sampledata, + const int* samplerow, + const int* samplecol, + const unsigned int* MDdata, + const int* MDrow, + const int* MDcol, + const int* burnin, + const int* mcmc, + const int* thin, + const int *uselecuyer, + const int *seedarray, + const int *lecuyerstream, + const int* verbose, + const double* thetastartdata, + const int* thetastartrow, + const int* thetastartcol, + const double* gammastartdata, + const int* gammastartrow, + const int* gammastartcol, + const double* tunevalue, + const double* thetaeqdata, + const int* thetaeqrow, + const int* thetaeqcol, + const double* thetaineqdata, + const int* thetaineqrow, + const int* thetaineqcol, + const int* storegamma, + const int* storetheta, + double* gammaacceptrate){ + + + // put together matrices + const Matrix MD(*MDrow, *MDcol, MDdata); + Matrix<> theta(*thetastartrow, *thetastartcol, thetastartdata); + Matrix<> gamma(*gammastartrow, *gammastartcol, gammastartdata); + const Matrix<> theta_eq(*thetaeqrow, *thetaeqcol, thetaeqdata); + const Matrix<> theta_ineq(*thetaineqrow, *thetaineqcol, thetaineqdata); + const int samplesize = (*samplerow) * (*samplecol); + + + MCMCPACK_PASSRNG2MODEL(MCMCpaircompare2d_impl, MD, theta, gamma, + theta_eq, theta_ineq, *tunevalue, + *burnin, *mcmc, *thin, + *verbose, *storegamma, *storetheta, + sampledata, samplesize, gammaacceptrate); + } +} + + +#endif + + + + + + + + + + diff -Nru r-cran-mcmcpack-1.5-0/src/cMCMCpaircompare2dDP.cc r-cran-mcmcpack-1.6-0/src/cMCMCpaircompare2dDP.cc --- r-cran-mcmcpack-1.5-0/src/cMCMCpaircompare2dDP.cc 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/src/cMCMCpaircompare2dDP.cc 2021-10-06 01:28:24.000000000 +0000 @@ -0,0 +1,936 @@ +////////////////////////////////////////////////////////////////////////// +// MCMCpaircompare.cc is C++ code to estimate a pairwise comparison model. +// +// KQ 3/18/2015 +// +// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn +// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn, +// and Jong Hee Park +////////////////////////////////////////////////////////////////////////// + +#ifndef MCMCPAIRCOMPARE2DDP_CC +#define MCMCPAIRCOMPARE2DDP_CC + +#include + +#include "MCMCrng.h" +#include "MCMCfcds.h" +#include "matrix.h" +#include "distributions.h" +#include "stat.h" +#include "la.h" +#include "ide.h" +#include "smath.h" +#include "rng.h" +#include "mersenne.h" +#include "lecuyer.h" + +#include // needed to use Rprintf() +#include // needed to allow user interrupts + +using namespace std; +using namespace scythe; + +/* MCMCpaircompare2dDP implementation. */ + +// update latent Ystar values in 2-d paired comparisons model +template +void paircompare2dDP_Ystar_update (Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& theta, const Matrix<>& gamma, + rng& stream){ + + const unsigned int N = MD.rows(); + for (unsigned int i = 0; i < N; ++i){ + const double gamma_i = gamma(MD(i,0)); + const double mean = std::cos(gamma_i)*theta(MD(i, 1), 0) + std::sin(gamma_i)*theta(MD(i, 1), 1) - + std::cos(gamma_i)*theta(MD(i, 2), 0) - std::sin(gamma_i)*theta(MD(i, 2), 1); + + const bool above = MD(i,1) == MD(i,3); // cand 1 chosen + const bool below = MD(i,2) == MD(i,3); // cand 2 chosen + if (above){ + Ystar(i) = stream.rtbnorm_combo(mean, 1.0, 0.0); + } + else if (below){ + Ystar(i) = stream.rtanorm_combo(mean, 1.0, 0.0); + } + else{ + Ystar(i) = stream.rnorm(mean, 1.0); + } + } +} + + + + + +// update theta values in 2-d paired comparisons model +template +void paircompare2dDP_theta_update (Matrix<>& theta, const Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& gamma, + const Matrix& theta_n, + const Matrix<>& theta_eq, + const Matrix<>& theta_ineq, + const vector< vector < double* > >& theta_Ystar_ptr, + const vector< vector < double* > >& theta_gamma_ptr, + const vector< vector < vector < double* > > >& theta_theta_ptr, + const vector< vector < double > >& theta_sign, + rng& stream){ + + const unsigned int J = theta.rows(); + const Matrix<> I = eye(2); + for (unsigned int j = 0; j < J; ++j){ + Matrix<> X(theta_n[j],2); + Matrix<> z(theta_n[j],1); + + // case 1: no equality constraints at all + if (theta_eq(j,0) == -999 && theta_eq(j,1) == -999){ + for (unsigned int i = 0; i < theta_n[j]; ++i){ + X(i,0)= std::cos(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + X(i,1)= std::sin(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + z(i,0) = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + (std::cos(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][0]+ + std::sin(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][1]); + } + Matrix<> v = invpd(crossprod(X) + I); + Matrix<> m = v * t(X)*z; + Matrix<> v_C = cholesky(v); + // case 1a: no inequality constraints at all + if (theta_ineq(j,0) == 0 && theta_ineq(j,1) == 0){ + Matrix<> theta_j_free = gaxpy(v_C, stream.rnorm(2, 1, 0, 1),m); + theta(j, _) = t(theta_j_free); + } + else{ // case 1b: some inequality constraint holds + Matrix<> theta_j_free = gaxpy(v_C, stream.rnorm(2, 1, 0, 1),m); + double ineq_holds = 0.0; + for (unsigned int jj = 0; jj < 2; ++jj){ + ineq_holds = std::min(ineq_holds, + theta_ineq(j,jj) * theta_j_free(jj)); + } + while (ineq_holds < 0){ + theta_j_free = gaxpy(v_C, stream.rnorm(2, 1, 0, 1), m); + double test = 0; + for (unsigned int jj = 0; jj < 2; ++jj) { + test = std::min(test, theta_ineq(j,jj) * theta_j_free(jj)); + } + ineq_holds = test; + } + theta(j,_) = t(theta_j_free); + } + } + else if (theta_eq(j,0) != -999 && theta_eq(j,1) == -999){ + // case 2: equality constraint on theta_j1 but not on theta_j2 + double xx = 0.0; + double xz = 0.0; + for (unsigned int i = 0; i < theta_n[j]; ++i){ + const double x_ji = std::sin(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + const double z_ji = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + (std::cos(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][0]+ + std::sin(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][1]) + -1*theta_sign[j][i] * theta_eq(j,0) * std::cos(*theta_gamma_ptr[j][i]); + xx += x_ji * x_ji; + xz += x_ji * z_ji; + } + const double v = 1.0/(xx + 1.0); + const double m = v * (xz); + // case 2a: no inequality constraint on theta_j2 + if (theta_ineq(j,1) == 0){ + theta(j,1) = stream.rnorm(m, std::sqrt(v)); + } + else{ // case 2b: inequality constraint on theta_j2 + if (theta_ineq(j,1) > 0) { // theta_j2 > 0 + theta(j,1) = stream.rtbnorm_combo(m, v, 0); + } else { // theta_j2 < 0 + theta(j,1) = stream.rtanorm_combo(m, v, 0); + } + } + } + else if (theta_eq(j,0) == -999 && theta_eq(j,1) != -999){ + // case 3: equality constraint on theta_j2 but not on theta_j1 + double xx = 0.0; + double xz = 0.0; + for (unsigned int i = 0; i < theta_n[j]; ++i){ + const double x_ji = std::cos(*theta_gamma_ptr[j][i]) * theta_sign[j][i]; + const double z_ji = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + (std::cos(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][0]+ + std::sin(*theta_gamma_ptr[j][i]) * *theta_theta_ptr[j][i][1]) + -1*theta_sign[j][i] * theta_eq(j,1) * std::sin(*theta_gamma_ptr[j][i]); + xx += x_ji * x_ji; + xz += x_ji * z_ji; + } + const double v = 1.0/(xx + 1.0); + const double m = v * (xz); + // case 3a: no inequality constraint on theta_j1 + if (theta_ineq(j,0) == 0){ + theta(j,0) = stream.rnorm(m, std::sqrt(v)); + } + else{ // case 3b: inequality constraint on theta_j1 + if (theta_ineq(j,0) > 0) { // theta_j1 > 0 + theta(j,0) = stream.rtbnorm_combo(m, v, 0); + } else { // theta_j1 < 0 + theta(j,0) = stream.rtanorm_combo(m, v, 0); + } + } + } + else if (theta_eq(j,0) != -999 && theta_eq(j,1) != -999){ + // case 4: equality constraints on both theta_j1 and theta_j2 + theta(j,0) = theta_eq(j,0); + theta(j,1) = theta_eq(j,1); + } + + } // end j loop + +} + +// update cluster_gamma in 2-d Dirichlet Process paired comparisons model +template +void paircompare2dDP_cluster_gamma_update (//Matrix<>& gamma, + const Matrix& gamma_n, + const vector< vector < double* > >& gamma_Ystar_ptr, + const vector< vector < vector< double* > > >& gamma_theta1_ptr, + const vector< vector < vector< double* > > >& gamma_theta2_ptr, + const double& tune, // A random walk step in the MH sampler : step \sim runif(-tune, tune). + const unsigned int& clustermcmc, + vector& gamma_trial, + vector& gamma_accept, + vector& judge_cluster_membership, //a vector of length I that stores judge i's cluster_membership + vector& cluster_gamma,//a vector of length K that stores the gamma value that charaterizes cluster k + vector& cluster_size, + rng& stream){ + + const unsigned int I = judge_cluster_membership.size(); + const unsigned int K = cluster_gamma.size(); + for (unsigned int k = 0; k < K; ++k){ + // if cluster k has no member for now, we generate the new k from its prior + if(cluster_size[k]==0){ + //I use 1.5707959999999 to denote PI/2 + cluster_gamma[k]=stream.runif()*1.5707959999999;//directly sampling from unif(0,PI/2) + }else{ //if cluster k has current members, we generate the new k through HM sampler + // we update gamma_k with the last value of this smaller HM sampler. + double gamma_k=cluster_gamma[k]; + for(unsigned int t=0;t1.5707959999999){ + gamma_k_new=gamma_k+(-2*stream.runif()+1)*tune; + } + //I assume a uniform prior for gamma, so I only consider likelihood for MH rejection rate. + double log_lik_old=0; + double log_lik_new=0; + double eta_i_j=0;//eta_i_j is the linear function value within the Probit link for judge i's j'th comparison. + double eta_i_j_new=0; + for(unsigned int i=0;i +void paircompare2dDP_cluster_weight_log_update ( + //vector& cluster_gamma,//a vector of length K that stores the gamma value that charaterizes cluster k + vector& cluster_weight_log,//on log scale + vector& cluster_size, + const double& alpha, + const unsigned int& I, + rng& stream){ + //IMPORTANT: cluster_weight_log is on the log scale + + + //work on cluster 1 through K-1: + const unsigned int K=cluster_weight_log.size(); + double cluster_cumulative_size= (double) I; + double weight_log=0; + double V_log=0; + double cumulative_one_minue_V_log=0;//must be initialized at zero + double rbeta_param1=0; + double rbeta_param2=0; + double beta_rv=0; + for(unsigned int k=0;k0.9999){ + beta_rv=stream.rbeta(rbeta_param1,rbeta_param2); + } + V_log=std::log(beta_rv); + weight_log=V_log+cumulative_one_minue_V_log; + cumulative_one_minue_V_log+=std::log(1-beta_rv); + cluster_weight_log[k]=weight_log; + } + //work on cluster K indexed by K-1 + cluster_weight_log[K-1]=cumulative_one_minue_V_log; +} + + + +// x: n array of original indices running from 0 to (n-1) +// n: length of x +// sample one element from x with equal likelihood +template +unsigned int SampleOneEqualLikelihood(int n, vector& x, rng& stream) { + //for (int i = 0; i < n; i++) + //x[i] = i; + //for (int i = 0; i < k; i++) { + double u = stream.runif();//Rcpp::runif(1, 0, 1)[0];//runif(1,0,1) has a NumericVecor output, so I need to index to get the double number out + int j = static_cast(n * u); + return x[j]; + // y[i] = x[j]; + // x[j] = x[--n]; + //} +} + +// draws 1 element from a const vector x with prob weights +// in array p +// n is length of x and p +template +int ProbSamp(const vector& x, const vector< double >& p, + const unsigned int& n, rng& stream) { + const double u = stream.runif(); + double cumprob = p[0]; + for (unsigned int i = 0; i < (n - 1); ++i) { + if (u <= cumprob) { + return x[i]; + } + else { + cumprob += p[i + 1]; + } + } + return x[n - 1]; +} + + +// update judge_cluster_membership in 2-d Dirichlet Process paired comparisons model +template +void paircompare2dDP_judge_cluster_membership_update (//Matrix<>& gamma, + const Matrix& gamma_n, + const vector< vector < double* > >& gamma_Ystar_ptr, + const vector< vector < vector< double* > > >& gamma_theta1_ptr, + const vector< vector < vector< double* > > >& gamma_theta2_ptr, + vector& judge_cluster_membership, //a vector of length I that stores judge i's cluster_membership + vector& cluster_gamma,//a vector of length K that stores the gamma value that charaterizes cluster k + vector& cluster_weight_log,//on log scale + vector& cluster_size, + const vector& cluster_labels, + unsigned int& unique_cluster, + rng& stream){ + + //IMPORTANT: I use the log-sum-exp trick to control underflow or overflow + // I store log-scale cluster probabiltiy values in cluster_probability_log. + // Moreover, I substract the max log-probabiltiy from every log-probability + + const unsigned int I = judge_cluster_membership.size(); + const unsigned int K = cluster_weight_log.size(); + vector< double > cluster_probability; + vector< double > cluster_probability_log; + cluster_probability.reserve(K); + cluster_probability_log.reserve(K); + + + + + for (unsigned int i = 0; i < I; ++i){ + double max_cluster_prob_log=-100000000000; + unsigned int cluster_i_new=0; + double denominator=0; + for(unsigned int k=0; k< K; ++k){ + //I assume a uniform prior for gamma, so I only consider likelihood for MH rejection rate. + double log_lik_i_k=0; //the log_lik of judge i being in cluster k + double eta_i_j=0;//eta_i_j is the linear function value within the Probit link for judge i's j'th comparison. + double gamma_i=cluster_gamma[k];//we let gamma_i try on each value in cluster_gamma + double cluster_prob_log=0; + for (unsigned int j = 0; j < gamma_n[i]; ++j){ + eta_i_j=std::cos(gamma_i)* *gamma_theta1_ptr[i][j][0] + std::sin(gamma_i)**gamma_theta1_ptr[i][j][1] - + std::cos(gamma_i)* *gamma_theta2_ptr[i][j][0] - std::sin(gamma_i)* *gamma_theta2_ptr[i][j][1] ; + log_lik_i_k += lndnorm(*gamma_Ystar_ptr[i][j], eta_i_j, 1);//lndnorm is the log scale normal density + + } + cluster_prob_log=cluster_weight_log[k]+log_lik_i_k;//both on log scale + cluster_probability_log[k]=cluster_prob_log; + if(cluster_prob_log>max_cluster_prob_log){ + max_cluster_prob_log=cluster_prob_log; + } + }//finish loop for computing cluster_prob_log for judge i + + //do the loop below to prevent underflow or overflow + for(unsigned int kk=0; kk< K; ++kk){ + cluster_probability_log[kk]-=max_cluster_prob_log; + } + //compute the cluster probability + + for(unsigned int kkk=0; kkk< K; ++kkk){ + denominator+=std::exp(cluster_probability_log[kkk]); + } + + + for(unsigned int kkkk=0; kkkk< K; ++kkkk){ + cluster_probability[kkkk]=std::exp(cluster_probability_log[kkkk])/denominator; + } + + + + cluster_i_new = ProbSamp(cluster_labels, cluster_probability, K, stream); + if(cluster_i_new != judge_cluster_membership[i]){ + + //bookkeeping for the number of unique clusters: + if(cluster_size[cluster_i_new]==0 && cluster_size[judge_cluster_membership[i]]>1){ + //new cluster used to be empty, and the current cluster has more than one member + ++unique_cluster; + }else if(cluster_size[cluster_i_new]>0 && cluster_size[judge_cluster_membership[i]]==1){ + //new cluster already has members, and the current cluster will be empty + --unique_cluster; + }//else: new cluster already has members,and the current cluster has more than one member. + //Don't need to change unique_cluster + + cluster_size[judge_cluster_membership[i]]-=1; + cluster_size[cluster_i_new]+=1; + judge_cluster_membership[i]=cluster_i_new; + }//else: the new cluster label for i is the same as the old one, we don't need to do anything + } + +} + + + + +// update alpha in 2-d Dirichlet Process paired comparisons model +template +void paircompare2dDP_alpha_update ( + vector& cluster_weight_log,//on log scale. I use the last element to update alpha + double& alpha, + const unsigned int& K, + const double& a, + const double& b, + rng& stream){ + + // shape=a+K-1; + // rate=b-cluster_weight_log[K-1]; + // cluster_weight_log is already on the log scale + // Rprintf("\nshape is: %7.3f", a+K-1, " "); + // Rprintf("\nrate is: %7.3f", b-cluster_weight_log[K-1], " "); + alpha=stream.rgamma(a+K-1, b-cluster_weight_log[K-1]); + // if(alpha < 0.1){ + // alpha=stream.rgamma(a+K-1, b-cluster_weight_log[K-1]); + // } + +} + + + +template +void MCMCpaircompare2dDP_impl (rng& stream, + const Matrix& MD, + Matrix<>& theta, Matrix<>& gamma, + Matrix<>& cluster_gamma_mat, + Matrix& judge_cluster_membership_mat, + const Matrix<>& theta_eq, + const Matrix<>& theta_ineq, + const double tune, + const unsigned int burnin, + const unsigned int mcmc, + const unsigned int clustermcmc, + const unsigned int thin, + const unsigned int verbose, + const bool storegamma, const bool storetheta, + double* sampledata, + const unsigned int samplesize, + double* gammaacceptrate, + double alpha, + const unsigned int clustermax, + const int alphafixed, + const double a, const double b){ + + + // constants + const unsigned int N = MD.rows();//data matrix + const unsigned int J = theta.rows(); + const unsigned int I = gamma.rows(); + const unsigned int tot_iter = burnin + mcmc; + const unsigned int nsamp = mcmc / thin; + const unsigned int K = clustermax;//the maximum number of clusters + + // starting values for Ystar + Matrix<> Ystar(N,1); + + + // pre-compute what can be pre-computed + Matrix theta_n(J,1,true,0); // vector of 0s. Item j's total instances of being compared. + Matrix gamma_n(I,1,true,0); // vector of 0s. Judge i's total instances of making comparisons. + for (unsigned int i = 0; i < N; ++i){ + gamma_n(MD(i,0)) += 1; + theta_n(MD(i,1)) += 1; + theta_n(MD(i,2)) += 1; + } + + vector< vector< double* > > theta_Ystar_ptr; + vector< vector< double* > > theta_gamma_ptr; + vector< vector< vector< double* > > > theta_theta_ptr;//the inner vector has two pointers for theta1 and theta2 + vector< vector< double > > theta_sign; + + vector< vector< double* > > gamma_Ystar_ptr; + vector< vector< vector< double* > > > gamma_theta1_ptr; //the inner vector has two pointers for theta1 and theta2 + vector< vector< vector< double* > > > gamma_theta2_ptr; //the inner vector has two pointers for theta1 and theta2 + vector judge_cluster_membership;//a vector of length I that stores judge i's cluster_membership + vector cluster_gamma;//a vector of length K that stores the gamma value that charaterizes cluster k + vector cluster_weight_log; + vector cluster_size; + vector cluster_labels; + + + //vector gamma_accept_rate_all; + vector gamma_trial; + vector gamma_accept; + + theta_Ystar_ptr.reserve(J); + theta_gamma_ptr.reserve(J); + theta_theta_ptr.reserve(J); + theta_sign.reserve(J); + + gamma_Ystar_ptr.reserve(I); + gamma_theta1_ptr.reserve(I); + gamma_theta2_ptr.reserve(I); + judge_cluster_membership.reserve(I); + cluster_gamma.reserve(K); + cluster_weight_log.reserve(K); + cluster_size.reserve(K); + + + //gamma_trial and gamma_accept record MH sampler proposal and accept for each cluster + gamma_trial.reserve(I); + gamma_accept.reserve(I); + + + for(unsigned int k=0; k( judge_cluster_membership.begin(), judge_cluster_membership.end() ).size(); + + + for(unsigned int k=0; k Ystar_j_ptr; + vector< double* > gamma_j_ptr; + vector< vector > theta_j_ptr; + vector< double > sign_j; + + Ystar_j_ptr.reserve(theta_n(j)); + gamma_j_ptr.reserve(theta_n(j)); + theta_j_ptr.reserve(theta_n(j)); + sign_j.reserve(theta_n(j)); + /* + Allocate proper length for the empty vectors (Ystar_j_ptr, gamma_j_ptr, theta_j_ptr, sign_j) + Then we save these pre-allocated vectors inside the outer layer vectors. + */ + theta_Ystar_ptr.push_back(Ystar_j_ptr); + theta_gamma_ptr.push_back(gamma_j_ptr); + theta_theta_ptr.push_back(theta_j_ptr); + theta_sign.push_back(sign_j); + } + + for (unsigned int i = 0; i < I; ++i){ + vector< double* > Ystar_i_ptr; + vector< vector< double* > > theta1_i_ptr; + vector< vector< double* > > theta2_i_ptr; + + Ystar_i_ptr.reserve(gamma_n(i)); + theta1_i_ptr.reserve(gamma_n(i)); + theta2_i_ptr.reserve(gamma_n(i)); + // choice_i.reserve(gamma_n(i)); + /* + Allocate proper length for the empty vectors (Ystar_j_ptr, theta1_i_ptr, theta2_i_ptr) + Then we save these pre-allocated vectors inside the outer layer vectors. + */ + gamma_Ystar_ptr.push_back(Ystar_i_ptr); + gamma_theta1_ptr.push_back(theta1_i_ptr); + gamma_theta2_ptr.push_back(theta2_i_ptr); + } + + + for (unsigned int i = 0; i < N; ++i){ + unsigned int resp = MD(i,0); + unsigned int c1 = MD(i,1); + unsigned int c2 = MD(i,2); + /* + Assign the starting values of Ystar, gamma, theta and sign to the above vectors. + This is convenient for accessing wanted objects in later updation. + */ + theta_Ystar_ptr[c1].push_back(&Ystar(i,0)); + theta_gamma_ptr[c1].push_back(&gamma(MD(i,0))); + vector theta_temp1; + theta_temp1.push_back(&theta(MD(i,2),0)); + theta_temp1.push_back(&theta(MD(i,2),1)); + theta_theta_ptr[c1].push_back(theta_temp1); + theta_sign[c1].push_back(1.0); + + theta_Ystar_ptr[c2].push_back(&Ystar(i,0)); + theta_gamma_ptr[c2].push_back(&gamma(MD(i,0))); + vector theta_temp2; + theta_temp2.push_back(&theta(MD(i,1),0)); + theta_temp2.push_back(&theta(MD(i,1),1)); + theta_theta_ptr[c2].push_back(theta_temp2); + theta_sign[c2].push_back(-1.0); + + gamma_Ystar_ptr[resp].push_back(&Ystar(i,0)); + gamma_theta1_ptr[resp].push_back(theta_temp2);//theta_temp2 represents c1 + gamma_theta2_ptr[resp].push_back(theta_temp1);//theta_temp1 represents c2 + + theta_temp1.clear(); + theta_temp2.clear(); + } + + + + + + // storage matrices (col major order) + Matrix<> theta_store; + Matrix<> gamma_store; + Matrix gamma_cluster_store; + Matrix unique_cluster_store; + Matrix<> alpha_store; + if (storetheta) + theta_store = Matrix<>(nsamp, J*2); + + if (storegamma){ + gamma_store = Matrix<>(nsamp, I); + gamma_cluster_store= Matrix(nsamp, I); + unique_cluster_store= Matrix(nsamp, 1); + } + + if(alphafixed==0){ + alpha_store= Matrix<>(nsamp, 1); + } + + + + + + + unsigned int count = 0; + // MCMC sampling occurs in this for loop + for (unsigned int iter = 0; iter < tot_iter; ++iter){ + + + // sample Ystar + paircompare2dDP_Ystar_update(Ystar, MD, theta, gamma, stream); + + + paircompare2dDP_cluster_gamma_update( + gamma_n, + gamma_Ystar_ptr, + gamma_theta1_ptr, + gamma_theta2_ptr, + tune, // A random walk step in the MH sampler : step \sim runif(-tune, tune). + clustermcmc, + gamma_trial, + gamma_accept, + judge_cluster_membership, //a vector of length I that stores judge i's cluster_membership + cluster_gamma,//a vector of length K that stores the gamma value that charaterizes cluster k + cluster_size, + stream); + + paircompare2dDP_cluster_weight_log_update ( + cluster_weight_log, + cluster_size, + alpha, + I, + stream); + + paircompare2dDP_judge_cluster_membership_update ( + gamma_n, + gamma_Ystar_ptr, + gamma_theta1_ptr, + gamma_theta2_ptr, + judge_cluster_membership, //a vector of length I that stores judge i's cluster_membership + cluster_gamma,//a vector of length K that stores the gamma value that charaterizes cluster k + cluster_weight_log,//on log scale + cluster_size, + cluster_labels, + unique_cluster, + stream); + + // update the gamma_i's according to the new judge_cluster_membership + for(unsigned int i=0;i candidates; + // for(unsigned int k=0;k 0 && iter % verbose == 0) { + Rprintf("\n\nMCMCpaircompare2dDP iteration %i of %i \n", (iter+1), tot_iter); + double gamma_accept_min = 1.0; + double gamma_accept_max = 0.0; + double gamma_accept_num = 0.0; + double gamma_accept_denom = 0.0; + for (unsigned int i=0; i gamma_accept_max) ? local_accrate : gamma_accept_max; + gamma_accept_num += local_accrate; + gamma_accept_denom += 1.0; + } + Rprintf("Minimum gamma acceptance rate: %7.3f", gamma_accept_min, " "); + Rprintf("\nMean gamma acceptance rate: %7.3f", gamma_accept_num/gamma_accept_denom, " "); + Rprintf("\nMaximum gamma acceptance rate: %7.3f", gamma_accept_max, " "); + Rprintf("\n\n"); + + } + + // store results + if ((iter >= burnin) && ((iter % thin == 0))) { + + // store theta + if (storetheta) { + for (unsigned int j = 0; j < J; ++j) { + theta_store(count, j) = theta(j,0); + theta_store(count, J+j) = theta(j, 1); + } + } + + + // store gamma + if (storegamma){ + gamma_store(count, _) = gamma; + for(unsigned int i=0; i output; + if(storetheta && ! storegamma) { + if(alphafixed==0){ + output = cbind(theta_store, alpha_store); + }else{ + output = theta_store; + } + + } else if (storegamma && ! storetheta){ + Matrix<> output_temp1; + Matrix<> output_temp2; + output_temp1 = cbind(gamma_cluster_store, unique_cluster_store); + output_temp2 = cbind(gamma_store, output_temp1); + + if(alphafixed==1){ + output = output_temp2; + }else{ + output = cbind(output_temp2, alpha_store); + } + + } else { + Matrix<> output_temp1; + Matrix<> output_temp2; + Matrix<> output_temp3; + output_temp1 = cbind(theta_store, gamma_store); + output_temp2 = cbind(gamma_cluster_store, unique_cluster_store); + output_temp3 = cbind(output_temp1, output_temp2); + if(alphafixed==1){ + output = output_temp3;//cbind(gamma_store, alpha_store); + }else{ + output = cbind(output_temp3, alpha_store); + } + + } + + for (unsigned int i = 0; i < samplesize; ++i) + sampledata[i] = output[i]; + + for (unsigned int i = 0; i < I; ++i) + gammaacceptrate[i] = gamma_accept[i]/gamma_trial[i]; + +} // end MCMCpaircompare implementation + + + + + + + + + +extern "C" { + + void + cMCMCpaircompare2dDP(double* sampledata, + const int* samplerow, + const int* samplecol, + const unsigned int* MDdata, + const int* MDrow, + const int* MDcol, + const int* burnin, + const int* mcmc, + const int* clustermcmc, + const int* thin, + const int *uselecuyer, + const int *seedarray, + const int *lecuyerstream, + const int* verbose, + const double* thetastartdata, + const int* thetastartrow, + const int* thetastartcol, + const double* gammastartdata, + const int* gammastartrow, + const int* gammastartcol, + const double* clustergammastartdata, + const int* clustergammastartrow, + const int* clustergammastartcol, + const int* judgeclustermembershipstartdata, + const int* judgeclustermembershipstartrow, + const int* judgeclustermembershipstartcol, + const double* tunevalue, + const double* thetaeqdata, + const int* thetaeqrow, + const int* thetaeqcol, + const double* thetaineqdata, + const int* thetaineqrow, + const int* thetaineqcol, + const int* storegamma, + const int* storetheta, + double* gammaacceptrate, + double* alpha, + const unsigned int* clustermax, + const int* alphafixed, + const double* a, + const double* b){ + + + // put together matrices + const Matrix MD(*MDrow, *MDcol, MDdata); + Matrix<> theta(*thetastartrow, *thetastartcol, thetastartdata); + Matrix<> gamma(*gammastartrow, *gammastartcol, gammastartdata); + Matrix<> cluster_gamma_mat(*clustergammastartrow, *clustergammastartcol, clustergammastartdata); + Matrix judge_cluster_membership_mat(*judgeclustermembershipstartrow, *judgeclustermembershipstartcol, judgeclustermembershipstartdata); + const Matrix<> theta_eq(*thetaeqrow, *thetaeqcol, thetaeqdata); + const Matrix<> theta_ineq(*thetaineqrow, *thetaineqcol, thetaineqdata); + const int samplesize = (*samplerow) * (*samplecol); + + + MCMCPACK_PASSRNG2MODEL(MCMCpaircompare2dDP_impl, MD, theta, gamma, cluster_gamma_mat, + judge_cluster_membership_mat, theta_eq, theta_ineq, *tunevalue, + *burnin, *mcmc, *clustermcmc, *thin, + *verbose, *storegamma, *storetheta, + sampledata, samplesize, gammaacceptrate, *alpha, + *clustermax, *alphafixed, *a, *b); + } +} + + +#endif + + + + + + + + + + diff -Nru r-cran-mcmcpack-1.5-0/src/cMCMCpaircompare.cc r-cran-mcmcpack-1.6-0/src/cMCMCpaircompare.cc --- r-cran-mcmcpack-1.5-0/src/cMCMCpaircompare.cc 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/src/cMCMCpaircompare.cc 2021-07-25 01:53:29.000000000 +0000 @@ -0,0 +1,315 @@ +////////////////////////////////////////////////////////////////////////// +// MCMCpaircompare.cc is C++ code to estimate a pairwise comparison model. +// +// KQ 3/18/2015 +// +// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn +// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn, +// and Jong Hee Park +////////////////////////////////////////////////////////////////////////// + +#ifndef MCMCPAIRCOMPARE_CC +#define MCMCPAIRCOMPARE_CC + +#include + +#include "MCMCrng.h" +#include "MCMCfcds.h" +#include "matrix.h" +#include "distributions.h" +#include "stat.h" +#include "la.h" +#include "ide.h" +#include "smath.h" +#include "rng.h" +#include "mersenne.h" +#include "lecuyer.h" + +#include // needed to use Rprintf() +#include // needed to allow user interrupts + +using namespace std; +using namespace scythe; + +/* MCMCpaircompare implementation. */ +template +void MCMCpaircompare_impl (rng& stream, + const Matrix& MD, + Matrix<>& theta, + Matrix<>& alpha, + const Matrix<>& theta_eq, + const Matrix<>& theta_ineq, + const double a0, + const double A0, + const int alphafixed, + const unsigned int burnin, + const unsigned int mcmc, + const unsigned int thin, + const unsigned int verbose, + const bool storealpha, + const bool storetheta, + double* sampledata, + const unsigned int samplesize){ + + + // constants + const unsigned int N = MD.rows(); + const unsigned int J = theta.rows(); + const unsigned int I = alpha.rows(); + const unsigned int tot_iter = burnin + mcmc; + const unsigned int nsamp = mcmc / thin; + + + // starting values for Ystar + Matrix<> Ystar(N,1); + + + // pre-compute what can be pre-computed + const double A0a0 = A0 * a0; + Matrix theta_n(J,1,true,0); // vector of 0s. Item j's total instances of being compared. + Matrix alpha_n(I,1,true,0); // vector of 0s. Judge i's total instances of making comparisons. + for (unsigned int i = 0; i < N; ++i){ + alpha_n(MD(i,0)) += 1; + theta_n(MD(i,1)) += 1; + theta_n(MD(i,2)) += 1; + } + + vector< vector< double* > > theta_Ystar_ptr; + vector< vector< double* > > theta_alpha_ptr; + vector< vector< double* > > theta_theta_ptr; + vector< vector< double > > theta_sign; + + vector< vector< double* > > alpha_Ystar_ptr; + vector< vector< double* > > alpha_theta1_ptr; + vector< vector< double* > > alpha_theta2_ptr; + + theta_Ystar_ptr.reserve(J); + theta_alpha_ptr.reserve(J); + theta_theta_ptr.reserve(J); + theta_sign.reserve(J); + + alpha_Ystar_ptr.reserve(I); + alpha_theta1_ptr.reserve(I); + alpha_theta2_ptr.reserve(I); + + for (unsigned int j = 0; j < J; ++j){ + vector< double* > Ystar_j_ptr; + vector< double* > alpha_j_ptr; + vector< double* > theta_j_ptr; + vector< double > sign_j; + + Ystar_j_ptr.reserve(theta_n(j)); + alpha_j_ptr.reserve(theta_n(j)); + theta_j_ptr.reserve(theta_n(j)); + sign_j.reserve(theta_n(j)); + /* + Allocate proper length for the empty vectors (Ystar_j_ptr, alpha_j_ptr, theta_j_ptr, sign_j) + Then we save these pre-allocated vectors inside the outer layer vectors. + */ + theta_Ystar_ptr.push_back(Ystar_j_ptr); + theta_alpha_ptr.push_back(alpha_j_ptr); + theta_theta_ptr.push_back(theta_j_ptr); + theta_sign.push_back(sign_j); + } + + for (unsigned int i = 0; i < I; ++i){ + vector< double* > Ystar_i_ptr; + vector< double* > theta1_i_ptr; + vector< double* > theta2_i_ptr; + + Ystar_i_ptr.reserve(alpha_n(i)); + theta1_i_ptr.reserve(alpha_n(i)); + theta2_i_ptr.reserve(alpha_n(i)); + /* + Allocate proper length for the empty vectors (Ystar_j_ptr, theta1_i_ptr, theta2_i_ptr) + Then we save these pre-allocated vectors inside the outer layer vectors. + */ + alpha_Ystar_ptr.push_back(Ystar_i_ptr); + alpha_theta1_ptr.push_back(theta1_i_ptr); + alpha_theta2_ptr.push_back(theta2_i_ptr); + } + + + for (unsigned int i = 0; i < N; ++i){ + unsigned int resp = MD(i,0); + unsigned int c1 = MD(i,1); + unsigned int c2 = MD(i,2); + /* + Assign the starting values of Ystar, alpha, theta and sign to the above vectors. + This is convenient for accessving wanted objects in later updation. + */ + theta_Ystar_ptr[c1].push_back(&Ystar(i,0)); + theta_alpha_ptr[c1].push_back(&alpha(MD(i,0))); + theta_theta_ptr[c1].push_back(&theta(MD(i,2))); + theta_sign[c1].push_back(1.0); + + theta_Ystar_ptr[c2].push_back(&Ystar(i,0)); + theta_alpha_ptr[c2].push_back(&alpha(MD(i,0))); + theta_theta_ptr[c2].push_back(&theta(MD(i,1))); + theta_sign[c2].push_back(-1.0); + + alpha_Ystar_ptr[resp].push_back(&Ystar(i,0)); + alpha_theta1_ptr[resp].push_back(&theta(MD(i,1))); + alpha_theta2_ptr[resp].push_back(&theta(MD(i,2))); + } + + + + + + // storage matrices (col major order) + Matrix<> theta_store; + Matrix<> alpha_store; + if (storetheta) + theta_store = Matrix<>(nsamp, J); + + if (storealpha) + alpha_store = Matrix<>(nsamp, I); + + + + unsigned int count = 0; + // MCMC sampling occurs in this for loop + for (unsigned int iter = 0; iter < tot_iter; ++iter){ + + // The following three functions are defined in MCMCfcds.h + // sample Ystar + paircompare_Ystar_update(Ystar, MD, theta, alpha, stream); + + // sample theta + paircompare_theta_update(theta, Ystar, MD, alpha, theta_n, theta_eq, + theta_ineq, + theta_Ystar_ptr, + theta_alpha_ptr, + theta_theta_ptr, + theta_sign, + stream); + + if (alphafixed == 0){ + // sample alpha + paircompare_alpha_update(alpha, Ystar, MD, theta, A0, A0a0, + alpha_n, + alpha_Ystar_ptr, + alpha_theta1_ptr, + alpha_theta2_ptr, + stream); + } + + // print results to screen + if (verbose > 0 && iter % verbose == 0) { + Rprintf("\n\nMCMCpaircompare iteration %i of %i \n", (iter+1), tot_iter); + //Rprintf("theta = \n"); + //for (int j=0; j= burnin) && ((iter % thin == 0))) { + + // store theta + if (storetheta) + theta_store(count, _) = theta; + + // store alpha + if (storealpha) + alpha_store(count, _) = alpha; + + count++; + } + + R_CheckUserInterrupt(); // allow user interrupts + + } // end MCMC sampling loop + + + + + // put output back into sampledata + Matrix<> output; + if(storetheta && ! storealpha) { + output = theta_store; + } else if (storealpha && ! storetheta){ + output = alpha_store; + } else { + output = cbind(theta_store, alpha_store); + } + + for (unsigned int i = 0; i < samplesize; ++i) + sampledata[i] = output[i]; + + +} // end MCMCpaircompare implementation + + + + + + + + + +extern "C" { + + void + cMCMCpaircompare(double* sampledata, + const int* samplerow, + const int* samplecol, + const unsigned int* MDdata, + const int* MDrow, + const int* MDcol, + const int* alphafixed, + const int* burnin, + const int* mcmc, + const int* thin, + const int *uselecuyer, + const int *seedarray, + const int *lecuyerstream, + const int* verbose, + const double* thetastartdata, + const int* thetastartrow, + const int* thetastartcol, + const double* astartdata, + const int* astartrow, + const int* astartcol, + const double* a0, + const double* A0, + const double* thetaeqdata, + const int* thetaeqrow, + const int* thetaeqcol, + const double* thetaineqdata, + const int* thetaineqrow, + const int* thetaineqcol, + const int* storealpha, + const int* storetheta){ + + + // put together matrices + const Matrix MD(*MDrow, *MDcol, MDdata); + Matrix<> theta(*thetastartrow, *thetastartcol, thetastartdata); + Matrix<> alpha(*astartrow, *astartcol, astartdata); + const Matrix<> theta_eq(*thetaeqrow, *thetaeqcol, thetaeqdata); + const Matrix<> theta_ineq(*thetaineqrow, *thetaineqcol, thetaineqdata); + const int samplesize = (*samplerow) * (*samplecol); + + + MCMCPACK_PASSRNG2MODEL(MCMCpaircompare_impl, MD, theta, alpha, + theta_eq, theta_ineq, *a0, *A0, + *alphafixed, + *burnin, *mcmc, *thin, + *verbose, *storealpha, *storetheta, + sampledata, samplesize); + } +} + + +#endif + + + + + + + + + + diff -Nru r-cran-mcmcpack-1.5-0/src/MCMCfcds.h r-cran-mcmcpack-1.6-0/src/MCMCfcds.h --- r-cran-mcmcpack-1.5-0/src/MCMCfcds.h 2019-11-18 19:55:23.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/src/MCMCfcds.h 2021-07-25 01:53:29.000000000 +0000 @@ -543,6 +543,137 @@ + +// update latent Ystar values in paired comparisons model +template +void paircompare_Ystar_update (Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& theta, const Matrix<>& alpha, + rng& stream){ + + const unsigned int N = MD.rows(); + for (unsigned int i = 0; i < N; ++i){ + const double alpha_i = alpha(MD(i,0)); + const double theta_1 = theta(MD(i,1)); + const double theta_2 = theta(MD(i,2)); + const double mean = alpha_i * (theta_1 - theta_2); + const bool above = MD(i,1) == MD(i,3); // cand 1 chosen + const bool below = MD(i,2) == MD(i,3); // cand 2 chosen + if (above){ + Ystar(i) = stream.rtbnorm_combo(mean, 1.0, 0.0); + } + else if (below){ + Ystar(i) = stream.rtanorm_combo(mean, 1.0, 0.0); + } + else{ + Ystar(i) = stream.rnorm(mean, 1.0); + } + } +} + + + + + +// update theta values in paired comparisons model +template +void paircompare_theta_update (Matrix<>& theta, const Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& alpha, + const Matrix& theta_n, + const Matrix<>& theta_eq, + const Matrix<>& theta_ineq, + const vector< vector < double* > >& theta_Ystar_ptr, + const vector< vector < double* > >& theta_alpha_ptr, + const vector< vector < double* > >& theta_theta_ptr, + const vector< vector < double > >& theta_sign, + rng& stream){ + + const unsigned int J = theta.rows(); + for (unsigned int j = 0; j < J; ++j){ + double xx = 0.0; + double xz = 0.0; + for (unsigned int i = 0; i < theta_n[j]; ++i){ + const double x_ji = theta_sign[j][i] * *theta_alpha_ptr[j][i]; + const double z_ji = *theta_Ystar_ptr[j][i] + theta_sign[j][i] * + *theta_alpha_ptr[j][i] * *theta_theta_ptr[j][i]; + /* + The reason for the above plus sign: -(-theta_sign[j][i])=+theta_sign[j][i]. + The other theta's sign in a comparison is the opposite to the theta in scrutiny. + */ + xx += x_ji * x_ji; + xz += x_ji * z_ji; + } + // prior for theta is N(0,1), therefore the posterior mean and variance are as follow + + const double v = 1.0/(xx + 1.0); + const double m = v * (xz); + // no equality constraints + if (theta_eq(j) == -999) { + if (theta_ineq(j) == 0) { // no inequality constraint + theta(j) = stream.rnorm(m, std::sqrt(v)); + } else if (theta_ineq(j) > 0) { // theta[j] > 0 + theta(j) = stream.rtbnorm_combo(m, v, 0); + } else { // theta[j] < 0 + theta(j) = stream.rtanorm_combo(m, v, 0); + } + } else { // equality constraints + theta(j) = theta_eq(j); + } + } + +} + + + + + + + + + +// update alpha values in paired comparisons model +template +void paircompare_alpha_update (Matrix<>& alpha, + const Matrix<>& Ystar, + const Matrix& MD, + const Matrix<>& theta, + const double& A0, + const double& A0a0, + const Matrix& alpha_n, + const vector< vector < double* > >& alpha_Ystar_ptr, + const vector< vector < double* > >& alpha_theta1_ptr, + const vector< vector < double* > >& alpha_theta2_ptr, + rng& stream){ + + const unsigned int I = alpha.rows(); + for (unsigned int i = 0; i < I; ++i){ + double xx = 0.0; + double xz = 0.0; + for (unsigned int j = 0; j < alpha_n[i]; ++j){ + const double x_ji = *alpha_theta1_ptr[i][j] - *alpha_theta2_ptr[i][j]; + const double z_ji = *alpha_Ystar_ptr[i][j]; + xx += x_ji * x_ji; + xz += x_ji * z_ji; + } + // prior for alpha is N(0,1), therefore the posterior mean and variance are as follow. Same as Bayesian regression + const double v = 1.0/(xx + A0); + const double m = v * (xz + A0a0); + alpha(i) = stream.rnorm(m, std::sqrt(v)); + } +} + + + + + + + + + + + + // factor analysis model with normal mean 0, precision F0 prior on // factor scores // X follows a multivariate normal distribution diff -Nru r-cran-mcmcpack-1.5-0/src/MCMCpack_init.c r-cran-mcmcpack-1.6-0/src/MCMCpack_init.c --- r-cran-mcmcpack-1.5-0/src/MCMCpack_init.c 2019-11-18 19:55:23.000000000 +0000 +++ r-cran-mcmcpack-1.6-0/src/MCMCpack_init.c 2021-07-25 01:53:29.000000000 +0000 @@ -45,6 +45,12 @@ void cHDPHMMnegbin(double *betaout, double *Pout, double *psout, double *sout, double *nuout, double *rhoout, double *tau1out, double *tau2out, int *comp1out, int *comp2out, double *sr1out, double *sr2out, double *mr1out, double *mr2out, double *gammaout, double *akout, double *thetaout, double *rhosizes, const double *Ydata, const int *Yrow, const int *Ycol, const double *Xdata, const int *Xrow, const int *Xcol, const int *K, const int *burnin, const int *mcmc, const int *thin, const int *verbose, const double *betastart, const double *Pstart, const double *nustart, const double *rhostart, const double *tau1start, const double *tau2start, const double *component1start, const double *gammastart, const double *akstart, const double *thetastart, const double *a_alpha, const double *b_alpha, const double *a_gamma, const double *b_gamma, const double *a_theta, const double *b_theta, const double *e, const double *f, const double *g, const double *rhostepdata, const int* uselecuyer, const int* seedarray, const int* lecuyerstream, const double *b0data, const double *B0data); extern void cHDPHMMpoisson(double *betaout, double *Pout, double *psout, double *sout, double *tau1out, double *tau2out, int *comp1out, int *comp2out, double *sr1out, double *sr2out, double *mr1out, double *mr2out, double *gammaout, double *akout, double *thetaout, const double *Ydata, const int *Yrow, const int *Ycol, const double *Xdata,const int *Xrow, const int *Xcol, const int *K, const int *burnin, const int *mcmc, const int *thin, const int *verbose, const double *betastart, const double *Pstart, const double *tau1start, const double *tau2start, const double *component1start, const double *gammastart, const double *akstart, const double *thetastart, const double *a_alpha, const double *b_alpha, const double *a_gamma, const double *b_gamma, const double *a_theta, const double *b_theta, const int* uselecuyer, const int* seedarray, const int* lecuyerstream, const double *b0data, const double *B0data); void cHDPHSMMnegbin(double *betaout, double *Pout, double *omegaout, double *sout, double *nuout, double *rhoout, double *tau1out, double *tau2out, int *comp1out, int *comp2out, double *sr1out, double *sr2out, double *mr1out, double *mr2out, double *gammaout, double *alphaout, double *rhosizes, const double *Ydata, const int *Yrow, const int *Ycol, const double *Xdata, const int *Xrow, const int *Xcol, const int *K, const int *burnin, const int *mcmc, const int *thin, const int *verbose, const double *betastart, const double *Pstart, const double *nustart, const double *rhostart, const double *tau1start, const double *tau2start, const double *component1start, const double *alphastart, const double *gammastart, const double *omegastart, const double *a_alpha, const double *b_alpha, const double *a_gamma, const double *b_gamma, const double *a_omega, const double *b_omega, const double *e, const double *f, const double *g, const double *r, const double *rhostepdata, const int* uselecuyer, const int* seedarray, const int* lecuyerstream, const double *b0data, const double *B0data); +void cMCMCpaircompare(double* sampledata, const int* samplerow, const int* samplecol, const unsigned int* MDdata, const int* MDrow, const int* MDcol, const int* alphafixed, const int* burnin, const int* mcmc, const int* thin, const int *uselecuyer, const int *seedarray, const int *lecuyerstream, const int* verbose, const double* thetastartdata, const int* thetastartrow, const int* thetastartcol, const double* astartdata, const int* astartrow, const int* astartcol, const double* a0, const double* A0, const double* thetaeqdata, const int* thetaeqrow, const int* thetaeqcol, const double* thetaineqdata, const int* thetaineqrow, const int* thetaineqcol, const int* storealpha, const int* storetheta); +void cMCMCpaircompare2d(double* sampledata, const int* samplerow, const int* samplecol, const unsigned int* MDdata, const int* MDrow, const int* MDcol, const int* burnin, const int* mcmc, const int* thin, const int *uselecuyer, const int *seedarray, const int *lecuyerstream, const int* verbose, const double* thetastartdata, const int* thetastartrow, const int* thetastartcol, const double* gammastartdata, const int* gammastartrow, const int* gammastartcol, const double* tunevalue, const double* thetaeqdata, const int* thetaeqrow, const int* thetaeqcol, const double* thetaineqdata, const int* thetaineqrow, const int* thetaineqcol, const int* storegamma, const int* storetheta, double* gammaacceptrate); +void cMCMCpaircompare2dDP(double* sampledata, const int* samplerow, const int* samplecol, const unsigned int* MDdata, const int* MDrow, const int* MDcol, const int* burnin, const int* mcmc, const int* clustermcmc, const int* thin, const int *uselecuyer, const int *seedarray, const int *lecuyerstream, const int* verbose, const double* thetastartdata, const int* thetastartrow, const int* thetastartcol, const double* gammastartdata, const int* gammastartrow, const int* gammastartcol, const double* clustergammastartdata, const int* clustergammastartrow, const int* clustergammastartcol, const int* judgeclustermembershipstartdata, const int* judgeclustermembershipstartrow, const int* judgeclustermembershipstartcol, const double* tunevalue, const double* thetaeqdata, const int* thetaeqrow, const int* thetaeqcol, const double* thetaineqdata, const int* thetaineqrow, const int* thetaineqcol, const int* storegamma, const int* storetheta, double* gammaacceptrate, double* alpha, const unsigned int* clustermax, const int* alphafixed, const double* a, const double* b); + + + /* .Call calls */ extern SEXP MCMClogituserprior_cc(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); @@ -92,6 +98,9 @@ {"cHDPHMMnegbin", (DL_FUNC) &cHDPHMMnegbin, 54}, {"cHDPHMMpoisson", (DL_FUNC) &cHDPHMMpoisson, 45}, {"cHDPHSMMnegbin", (DL_FUNC) &cHDPHSMMnegbin, 54}, + {"cMCMCpaircompare", (DL_FUNC) &cMCMCpaircompare, 30}, + {"cMCMCpaircompare2d", (DL_FUNC) &cMCMCpaircompare2d, 29}, + {"cMCMCpaircompare2dDP", (DL_FUNC) &cMCMCpaircompare2dDP, 41}, {NULL, NULL, 0} };