This function draws npv plausible values for each person from their posterior density.
PV(estobj, ...) # S3 method for fourpl PV(estobj, npv = 10, approx = TRUE, thinning = 6, burnin = 10, mult = 2, ...) # S3 method for gpcm PV(estobj, npv = 10, approx = TRUE, thinning = 6, burnin = 10, mult = 2, ...) # S3 method for gpcm4pl PV(estobj, npv = 10, approx = TRUE, thinning = 6, burnin = 10, mult = 2, ...) # S3 method for pv print(x, ...) # S3 method for pv summary(object, nrowmax = 15, ...)
| estobj | An object which originates from using |
|---|---|
| ... | More arguments |
| npv | The number of (effectively returned) plausible values - default is 10. |
| approx | Whether a normal approximation |
| thinning | 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). |
| burnin | 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. |
| mult | Multiplication constant (default = 2). Use this parameter to vary the width of the proposal distribution - which is |
| x | An object of class |
| object | An object of class |
| nrowmax | 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.
Manuel Reif
################# 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!#> PP Version: 0.6.3.11 #> #> Call: PV.fourpl(estobj = res2pleap) #> - job started @ Mon May 24 13:27:58 2021 #> #> 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.3.11 #> #> Call: PV.fourpl(estobj = res2pleap, approx = FALSE) #> - job started @ Mon May 24 13:27:58 2021 #> #> 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#> [1] 0.0334621629 -0.1829705082 0.0348347988 -0.1514444247 -0.0932831102 #> [6] -0.0421735214 -0.1183732851 0.0202321291 -0.0008060551 -0.0892993107#> [,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#> [,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#> [,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