This function draws npv plausible values for each person from their posterior density.
PV(estobj, ...)
# S3 method for class 'fourpl'
PV(estobj, npv = 10, approx = TRUE, thinning = 6, burnin = 10, mult = 2, ...)
# S3 method for class 'gpcm'
PV(estobj, npv = 10, approx = TRUE, thinning = 6, burnin = 10, mult = 2, ...)
# S3 method for class 'gpcm4pl'
PV(estobj, npv = 10, approx = TRUE, thinning = 6, burnin = 10, mult = 2, ...)
# S3 method for class 'pv'
print(x, ...)
# S3 method for class 'pv'
summary(object, nrowmax = 15, ...)An object which originates from using PP_4pl(), PP_gpcm() or PPall(). EAP estimation is strongly recommended (type = "eap"), when plausible values are drawn afterwards, because the EAP estimate is used as starting point for the MH algorithm.
More arguments
The number of (effectively returned) plausible values - default is 10.
Whether a normal approximation N(mu,sigma2) is used to draw the plausible values. Default = TRUE. If FALSE a Metropolitan-Hastings-Algorithm will draw the values.
A numeric vector of length = 1. If approx = FALSE, a Metropolitan-Hastings-Algorithm draws the plausible values. To avoid autocorrelation, thinning takes every kth value as effective plausible value. The default is 6 (every 6th value is taken), which works appropriately in almost all cases here (with default settings).
How many draws should be discarded at the chains beginning? Default is 10 - and this seems reasonable high (probably 5 will be enough as well), because starting point is the EAP.
Multiplication constant (default = 2). Use this parameter to vary the width of the proposal distribution - which is N(theta_v,mult*SE_eap) - when a MH-Algorithm is applied. So the constant quantifies the width in terms of multiples of the EAP standard error. 2 works fine with the default thinning. If the supplied value is large, thinning can take lower values without causing autocorrelation.
An object of class pv which is the result of using the PV() function
An object of class pv which is the result of using the PV() function
When printing the matrix of estimates - how many rows should be shown? Default = 15.
The function returns a list which main element is pvdraws. This is a matrix of size number_of_persons x npv - so if 10 plausible values are requested for 100 persons, a 100x10 matrix is returned.
Mislevy, R. J. (1991). Randomization-based inference about latent variables from complex samples. Psychometrika, 56(2), 177-196.
Von Davier, M., Gonzalez, E., & Mislevy, R. (2009). What are plausible values and why are they useful. IERI monograph series, 2, 9-36.
Kruschke, J. (2010). Doing Bayesian data analysis: A tutorial introduction with R. Academic Press.
################# Plausible values #############################################################
### 4 PL model ######
### data creation ##########
set.seed(1522)
# intercepts
diffpar <- seq(-3,3,length=12)
# slope parameters
sl <- round(runif(12,0.5,1.5),2)
la <- round(runif(12,0,0.25),2)
ua <- round(runif(12,0.8,1),2)
# response matrix
awm <- matrix(sample(0:1,10*12,replace=TRUE),ncol=12)
# EAP estimation - 2pl model
res2pleap <- PP_4pl(respm = awm,thres = diffpar, slopes = sl,type = "eap")
#> Warning: all mu's are set to 0!
#> Warning: all sigma2's are set to 1!
#> Estimating: 2pl model ...
#> type = eap
#> Estimation finished!
# draw 10 plausible values
res_pv <- PV(res2pleap)
summary(res_pv)
#> PP Version: 0.6.4.0
#>
#> Call: PV.fourpl(estobj = res2pleap)
#> - job started @ Tue Apr 21 13:42:32 2026
#>
#> Estimation type: Plausible values
#>
#> -------------------------------------
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0.6953 -0.4238 0.5376 0.1220 1.1666 -0.4869 0.3956 -0.1345 0.4868
#> [2,] 1.8993 1.3331 1.9698 2.6371 2.9233 1.6535 1.2706 2.0875 3.1040
#> [3,] -0.1233 0.3234 0.2542 1.1840 0.8564 1.2455 0.2902 2.0925 0.9092
#> [4,] -0.0955 0.4056 -0.0211 1.1267 -0.3397 1.2787 0.3276 0.0335 0.7451
#> [5,] -0.3687 0.4171 0.0624 0.9619 0.6330 1.1589 -0.1714 1.2052 1.0300
#> [6,] -0.1402 -0.8282 -0.0552 -0.6590 -0.2023 0.3184 -0.2779 -0.7942 -0.4154
#> [7,] 0.4481 -0.6536 1.5925 0.1980 0.4367 -0.5248 -0.5444 0.8554 -0.1203
#> [8,] -1.0657 -0.1015 -0.5556 -0.5727 -0.8175 0.2507 0.2714 -0.5924 -0.4320
#> [9,] -0.1300 0.8070 0.9653 1.0756 0.0884 -0.4184 0.7620 -0.3311 0.7093
#> [10,] 0.0168 -0.0284 -1.0015 -1.0178 -0.5916 -1.1864 -2.0544 -0.0423 0.1119
#> [,10]
#> [1,] 1.4409
#> [2,] 2.9403
#> [3,] 0.5839
#> [4,] -0.1149
#> [5,] 0.0241
#> [6,] 0.2265
#> [7,] 0.9147
#> [8,] -1.3009
#> [9,] -0.2065
#> [10,] -0.2729
# draw 10 plausible values - use a metropolitan hastings algorithm
res_pv2 <- PV(res2pleap,approx = FALSE)
summary(res_pv2)
#> PP Version: 0.6.4.0
#>
#> Call: PV.fourpl(estobj = res2pleap, approx = FALSE)
#> - job started @ Tue Apr 21 13:42:32 2026
#>
#> Estimation type: Plausible values
#>
#> -------------------------------------
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0.5176 0.6702 0.8384 0.1567 0.3942 0.3183 1.1856 0.4123 1.3934
#> [2,] 2.6354 1.7938 2.2411 0.9559 1.5905 1.8675 1.9793 2.1926 2.4647
#> [3,] 1.1811 0.9038 0.2391 0.3320 0.5533 0.0269 0.3355 0.7387 2.3213
#> [4,] 0.6674 1.2358 -0.1539 0.4765 -0.9007 0.0618 -0.0947 0.1080 0.3025
#> [5,] 1.2287 0.7294 1.1350 0.1083 0.2721 0.2760 0.1822 0.9192 -0.1466
#> [6,] -0.3054 0.4222 0.1953 -0.1967 -1.2530 -0.9176 0.8579 -1.0329 -0.1456
#> [7,] -0.3392 1.3096 0.8709 0.1195 1.0242 0.6261 1.1982 0.1842 0.1152
#> [8,] -1.2496 -1.1042 0.2793 0.6694 -1.1338 -0.2512 -0.3135 -0.2266 -0.8891
#> [9,] 0.9876 1.1833 -0.3425 0.2033 -0.2409 1.1824 -0.5005 0.4128 1.8582
#> [10,] 0.1786 -0.3631 -0.1823 -0.8970 0.0599 -1.4422 -0.2283 -1.0609 0.3647
#> [,10]
#> [1,] 0.9946
#> [2,] 2.4778
#> [3,] 0.8432
#> [4,] -0.5142
#> [5,] 1.0123
#> [6,] 0.7441
#> [7,] 0.7339
#> [8,] -0.7373
#> [9,] 0.4449
#> [10,] -1.4295
# ------ check the PVs
# -- autocorrelation?
autocor <- function(acv)
{
cor(acv[-1],acv[-length(acv)])
}
res_pvac <- PV(res2pleap,approx = FALSE,npv = 200)
# independent draws - so there cannot be any systematic autocorrelation when
# approx = TRUE. So this acts as a kind of benchmark for the MH-Alg.
res_pvac2 <- PV(res2pleap,approx = TRUE,npv = 200)
apply(res_pvac$pvdraws,1,autocor)
#> [1] -0.09218511 -0.09269937 -0.02118295 -0.10091570 -0.04938669 -0.04611125
#> [7] -0.04625685 -0.02323778 -0.10955832 -0.08582139
apply(res_pvac2$pvdraws,1,autocor)
#> [1] 0.0334621629 -0.1829705082 0.0348347988 -0.1514444247 -0.0932831102
#> [6] -0.0421735214 -0.1183732851 0.0202321291 -0.0008060551 -0.0892993107
# -- autocorrelation distr?
apply(res_pvac$pvdraws,1,quantile)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> 0% -1.1189126 0.6482584 -0.8599093 -1.17206073 -1.06657367 -1.9113879
#> 25% 0.1243905 1.6327885 0.2907127 -0.05658439 -0.03965872 -0.5590972
#> 50% 0.5966065 2.1110250 0.7791557 0.37603789 0.54105135 -0.1806396
#> 75% 1.1131146 2.6571856 1.1399531 0.73833154 1.01487874 0.1831070
#> 100% 2.7299487 3.7852102 2.1300067 2.24528325 1.92887105 1.7386137
#> [,7] [,8] [,9] [,10]
#> 0% -1.39552654 -2.9369528 -0.89186083 -2.5718982
#> 25% -0.03756925 -1.1410654 0.04793972 -1.2907359
#> 50% 0.44731563 -0.5981778 0.46023882 -0.9884915
#> 75% 0.93801834 -0.1137122 0.89671504 -0.4047435
#> 100% 2.07215601 0.9584107 2.37528014 0.9008073
apply(res_pvac2$pvdraws,1, quantile)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> 0% -1.91636101 0.01988437 -1.2642093 -1.29134152 -1.441746554 -1.8945768
#> 25% 0.06152027 1.66679080 0.2228406 -0.04999196 0.009426829 -0.5331227
#> 50% 0.65586226 2.08172297 0.6867953 0.36135697 0.552124816 -0.1308765
#> 75% 0.96952681 2.57572512 1.1293040 0.68278043 0.973946245 0.2960726
#> 100% 2.39979403 3.93395552 2.3552885 1.69208634 1.878495326 1.9703464
#> [,7] [,8] [,9] [,10]
#> 0% -1.40250925 -2.4268740 -1.38035907 -2.3991210
#> 25% 0.04597495 -0.9474394 0.02987099 -1.1518491
#> 50% 0.46163170 -0.5717532 0.46486795 -0.7977427
#> 75% 0.93354670 -0.1431859 0.88224389 -0.4279617
#> 100% 2.03832546 1.2638010 1.99662599 0.6953115
### GPCM model ######
# some threshold parameters
THRES <- matrix(c(-2,-1.23,1.11,3.48,1
,2,-1,-0.2,0.5,1.3,-0.8,1.5),nrow=2)
# slopes
sl <- c(0.5,1,1.5,1.1,1,0.98)
awmatrix <- matrix(c(1,0,2,0,1,1,1,0,0,1
,2,0,0,0,0,0,0,0,0,1,1,2,2,1,1,1,1,0,0,1),byrow=TRUE,nrow=5)
# EAP estimation
resgpcmeap <- PP_gpcm(respm = awmatrix,thres = THRES, slopes = sl,type = "eap")
#> Warning: all mu's are set to 0!
#> Warning: all sigma2's are set to 1!
#> Estimating: GPCM ...
#> type = eap
#> Estimation finished!
res_gpcmpv <- PV(resgpcmeap,approx = FALSE,npv = 20)
### GPCM and 4PL model ######
# some threshold parameters
THRES <- matrix(c(-2,-1.23,1.11,3.48,1
,2,-1,-0.2,0.5,1.3,-0.8,1.5),nrow=2)
# slopes
sl <- c(0.5,1,1.5,1.1,1,0.98)
THRESx <- THRES
THRESx[2,1:3] <- NA
# for the 4PL item the estimated parameters are submitted,
# for the GPCM items the lower asymptote = 0
# and the upper asymptote = 1.
la <- c(0.02,0.1,0,0,0,0)
ua <- c(0.97,0.91,1,1,1,1)
awmatrix <- matrix(c(1,0,1,0,1,1,1,0,0,1
,2,0,0,0,0,0,0,0,0,1
,1,2,2,1,1,1,1,0,0,1),byrow=TRUE,nrow=5)
model2est <- findmodel(THRESx)
# EAP estimation
respcmeap1 <- PPall(respm = awmatrix,thres = THRESx,
slopes = sl,lowerA = la, upperA=ua, type = "eap",
model2est=model2est)
#> Estimating: mixed 4PL, GPCM ...
#> type = eap
#> Estimation finished!
res_mixedpv_1 <- PV(respcmeap1,approx = FALSE,npv = 200)
# rowMeans of plausible values should approximate the EAPs
rowMeans(res_mixedpv_1$pvdraws)
#> [1] 0.093466412 -0.083368876 -1.607109734 1.298907846 0.004484088
# EAPs
respcmeap1
#> estimate SE
#> [1,] 0.0682 0.5905
#> [2,] -0.0653 0.5950
#> [3,] -1.5237 0.7046
#> [4,] 1.3112 0.6142
#> [5,] -0.0417 0.6119
# show the quantiles of the empirical distribution
apply(res_mixedpv_1$pvdraws,1,quantile)
#> [,1] [,2] [,3] [,4] [,5]
#> 0% -1.3485464 -1.57470138 -3.8011143 -0.01345653 -1.45492521
#> 25% -0.3082155 -0.54806041 -2.0808571 0.82804669 -0.50233232
#> 50% 0.1284492 -0.01708521 -1.5534572 1.27125852 0.00656254
#> 75% 0.4943450 0.33145978 -1.0425829 1.72900086 0.49265666
#> 100% 1.6374132 1.49828104 0.2286223 2.90438148 1.62675981