Introduction to conquestr

These functions allows users to include ‘ACER ConQuest’ within their R workflow. This is particularly useful for users of the Mac OS X version (and Windows console) version of ‘ACER ConQuest’ which does not support ploting.

This vignette demonstrates how to use the package built-in demonstration files to:

Call ConQuest and run a syntax file

First, ensure you have the most current version of ConQuest. Version > 5.17.4 is required. If you are unsure, check ConQuest on the ACER Shop. You will also need to install the conquestr library:

install.packages("conquestr")
library("conquestr")

You must pass the ConQuestCall function at least two pieces of information: the install location of ConQuest and a valid syntax file. conquestr has included syntax files that can be used to run a model and test other package functions. A default syntax file will be called is you do not explicitly provide one. conquestr will search for a valid installation of ‘ACER ConQuest’ unless the executable is explicitly declared. The search locations are (in order):

The default install location (Windows:%ProgramFiles%\ACER ConQuest\ConQuestConsole.exe and Mac OS X: /Applications/ConQuest/ConQuest).
The Program Files folder on Windows, and the Application folder on Mac OS X
The current user’s Desktop

# if you don't provide a syntax file, using the argumenmt `cqc=`, then the in-built demo syntax file will be run. 
ConQuestCall()

# the following output is produced:

ConQuest build: Mar 10 2021
<TYPE> Version
This version expires <DATE>
submit <PATH>conquestr/extdata/conquestabout.cqc;
=>dir;
<PATH>
=>about;
Developed by
    Australian Council for Educational Research
    University of California,Berkeley

Your key: <KEY>
Expires: <DATE>

<TYPE> Build:  <DATE>
Version: 5.17.4

Programmers
    Ray Adams,Margaret Wu,Greg Macaskill,Sam Haldane,Xiao Xun Sun,Dan Cloney
End of Program

When you use ConQuestCall to call ‘ACER ConQuest’ and run a syntax file, some options are set by default. The same options can be set within ‘ACER ConQuest’ by using the set command (set progress = yes, exit_on_error = yes, storecommands = yes, warnings = no; which is syntactically the same as using the helper function set conquestr = true;). WARNING - this will make it easy to overwrite output as you will not be prompted or warned. To turn off these settings your syntax file must explicitly change these using the set command. For more information see the set command.

When calling ‘ACER ConQuest’ from R, the ConQuest working directory will default to the current R working directory. This makes it easy to write portable syntax using relative paths. If you need ‘ACER ConQuest’ to use a different working directory, set it in your syntax file.

The next example runs a small analysis and creates a system file that is read into the R object myExSys.

# set up
oldWd <- getwd()
myTmpDir <- tempdir()
setwd(myTmpDir)
file.copy(from = c(
  system.file("extdata", "ex1.cqc", package = "conquestr"),
  system.file("extdata", "ex1.dat", package = "conquestr"),
  system.file("extdata", "mysysfile.cqs", package = "conquestr")
  ), to = myTmpDir)
#> [1] TRUE TRUE TRUE

# run ConQuest
# ConQuestCall is not run here, as it requires a local install of ConQuest 
# ConQuestCall(cqc = file.path(myTmpDir, "ex1.cqc"), stdout = NULL) # searches for valid install of CQ
myExSys <- ConQuestSys(myCqs = file.path(myTmpDir, "mysysfile.cqs"))
#> System file is compressed. Decompressing ...
#> complete

# see the content of the sysfile
str(myExSys)
#> List of 202
#>  $ compressedString       : chr "uncompressed"
#>  $ builddate              : chr "May 17 2025"
#>  $ writedate              : chr "Sun May 18 10:02:56 2025\n"
#>  $ cqs_version            : int 29
#>  $ gNCases                :List of 11
#>   ..$ : num 1000
#>   ..$ : num 1000
#>   ..$ : num 1000
#>   ..$ : num 1000
#>   ..$ : num 999
#>   ..$ : num 75
#>   ..$ : num 925
#>   ..$ : num 1000
#>   ..$ : num 999
#>   ..$ : num 925
#>   ..$ : num 1000
#>  $ gNDim                  : int 1
#>  $ gNGins                 : int 12
#>  $ gNPlausiblesEstimate   : int 5
#>  $ gMLEExist              : logi TRUE
#>  $ gWLEExist              : logi TRUE
#>  $ gEAPExist              : logi TRUE
#>  $ gPlausibleExist        : logi TRUE
#>  $ gSystemMissing         : num -1.8e+308
#>  $ gApplyFilter           : logi FALSE
#>  $ gFilter                :List of 4
#>   ..$ String : logi [1, 1] NA
#>   ..$ Items  : int 0
#>   ..$ Columns: int 0
#>   ..$ Size   : int 0
#>  $ gBeta                  : num [1, 1] 0
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gOldBeta               : num [1, 1] 0
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gBestBeta              : num [1, 1] 0
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gXsi                   : num [1:12, 1] -0.704 -1.251 -1.092 -0.231 0.111 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:12] "" "" "" "" ...
#>   .. ..$ : chr ""
#>  $ gOldXsi                : num [1:12, 1] -0.704 -1.251 -1.092 -0.231 0.111 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:12] "" "" "" "" ...
#>   .. ..$ : chr ""
#>  $ gBestXsi               : num [1:12, 1] -0.704 -1.251 -1.092 -0.231 0.111 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:12] "" "" "" "" ...
#>   .. ..$ : chr ""
#>  $ gTau                   : logi [1, 1] NA
#>  $ gOldTau                : logi [1, 1] NA
#>  $ gBestTau               : logi [1, 1] NA
#>  $ gQuickErrorsXsi        : num [1:12, 1] 0.00513 0.00605 0.00572 0.00479 0.00476 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr ""
#>  $ gQuickErrorsTau        : logi [1, 1] NA
#>  $ gQuickErrorsSigma      : num [1, 1] 0.0015
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gQuickErrorsBeta       : logi [1, 1] NA
#>  $ gMasterTheta           : num [1:2000, 1] 0.69 0.679 0.89 0.15 -0.626 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:2000] "" "" "" "" ...
#>   .. ..$ : chr ""
#>  $ gTheta                 : num [1:2000, 1] 0.455 0.445 0.641 -0.047 -0.769 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:2000] "" "" "" "" ...
#>   .. ..$ : chr ""
#>  $ gVariance              : num [1, 1] 0.865
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gOldVariance           : num [1, 1] 0.865
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gBestVariance          : num [1, 1] 0.865
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gHistoryWeights        : logi [1, 1] NA
#>  $ gOldHistory            : logi [1, 1] NA
#>  $ gBestHistory           : logi [1, 1] NA
#>  $ gYBetaAll              : num [1:1000, 1] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:1000] "" "" "" "" ...
#>   .. ..$ : chr ""
#>  $ gWtFactor              :List of 4
#>   ..$ : num 1
#>   ..$ : num 1
#>   ..$ : num 1
#>   ..$ : num 1
#>  $ gSuffXsi               : num [1:12, 1] -644 -743 -716 -548 -476 -764 -776 -766 -852 -776 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr ""
#>  $ gSuffTau               : logi [1, 1] NA
#>  $ gModelText             : chr "item"
#>  $ gFormatText            : chr " pid 1-5 responses 12-23 v1 12 v2 15"
#>  $ gRegressionText        : chr ""
#>  $ gGroupText             : chr "v1"
#>  $ gOSSCP                 : num [1, 1] 999
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr ""
#>   .. ..$ : chr ""
#>  $ gLOSSCP                : logi [1, 1] NA
#>  $ gLSSCP                 : logi [1, 1] NA
#>  $ gFullSSCP              : logi [1, 1] NA
#>  $ gFullSums              : logi [1, 1] NA
#>  $ gMinAlpha              : num 0
#>  $ gModelEstimated        : logi TRUE
#>  $ gIntegrationMethod     :'data.frame': 7 obs. of  2 variables:
#>   ..$ gIntegrationMethod    : int [1:7] 1 2 3 4 5 6 7
#>   ..$ gIntegrationMethodText: chr [1:7] "Bock Aitkin" "Monte Carlo" "Gauss-Hermite Quadrature" "Joint Maximum Likelihood" ...
#>  $ gPopulation            : int 1
#>  $ gSeeds                 : int 5
#>  $ gMaxSinceBests         : int 100
#>  $ gInnerLoops            : int 10
#>  $ gWarnings              : logi TRUE
#>  $ gEstsToLog             : logi TRUE
#>  $ gKeepLast              : logi FALSE
#>  $ gAddExtension          : logi TRUE
#>  $ gMLEMax                : num 15
#>  $ gPlotWinMax            : int 100
#>  $ gZero                  : num 0.3
#>  $ gRespMiss              : int 2
#>  $ gDatafileName          : chr "ex1.dat"
#>  $ gDatafileFormats       : int 0
#>  $ gDatafileNameDisplay   : chr "ex1.dat"
#>  $ gStopReason            : int 2
#>  $ gImplicit              :List of 2
#>   ..$ Name  :List of 1
#>   .. ..$ : chr "item"
#>   ..$ Levels:List of 1
#>   .. ..$ : int 12
#>  $ gNImpValue             : int 12
#>  $ gPIDVar                : int 0
#>  $ gModelVariables        :List of 2
#>   ..$ XV: list()
#>   ..$ IV:List of 1
#>   .. ..$ : int 0
#>  $ gNRec                  : int 1
#>  $ gResponseLookUp        :List of 2
#>   ..$ VarNumber: int -2
#>   ..$ Value    :List of 2
#>   .. ..$ : chr "0"
#>   .. ..$ : chr "1"
#>  $ gPreKeyLookUp          :List of 2
#>   ..$ VarNumber: int -2
#>   ..$ Value    :List of 6
#>   .. ..$ : chr "M"
#>   .. ..$ : chr "a"
#>   .. ..$ : chr "b"
#>   .. ..$ : chr "c"
#>   .. ..$ : chr "d"
#>   .. ..$ : chr "e"
#>  $ gNDataRecords          : int 11989
#>  $ gFacetVariables        :List of 2
#>   ..$ XV: list()
#>   ..$ IV: list()
#>  $ gGroupVariables        :List of 2
#>   ..$ XV:List of 1
#>   .. ..$ : int 1
#>   ..$ IV: list()
#>  $ gWeightVariable        :List of 2
#>   ..$ XV: list()
#>   ..$ IV: list()
#>  $ gTDFileV               :List of 2
#>   ..$ XV: list()
#>   ..$ IV:List of 1
#>   .. ..$ : int 0
#>  $ gValidC                : list()
#>  $ gFileRebuildNeeded     : logi FALSE
#>  $ gAMatrixImportFileName : chr ""
#>  $ gCMatrixImportFileName : chr ""
#>  $ gTitle                 : chr "ConQuest: Generalised Item Response Modelling Software"
#>  $ gSubmitMode            : logi TRUE
#>  $ gMaxCats               : int 2
#>  $ gConvergenceOK         : logi FALSE
#>  $ gParameterConvCriterion: num 1e-04
#>  $ gDevianceConvCriterion : num 1e-04
#>  $ gFitDraws              : int 5
#>  $ gMaxIterations         : int 1000
#>  $ gAccuracy              : int 15
#>  $ gPVNodes               : int 2000
#>  $ gFitNodes              : int 2000
#>  $ gIteration             : int 23
#>  $ gBestIter              : int 23
#>  $ gStdError              : int 0
#>  $ gIFit                  : logi TRUE
#>   [list output truncated]
#>  - attr(*, "class")= chr [1:2] "list" "conQuestSysFile"

# tidy up
file.remove(
  list.files(myTmpDir, all.files = TRUE, full.names = TRUE, pattern = "^myi|^myS|ex1"),
  recursive=TRUE
)
#> [1] TRUE TRUE
setwd(oldWd)

Read in a ConQuest system file

The above example created a system file. Alternatively, there is an inbuilt example of a system file.

conquestr can read in a system file, using the ConQuestSys function. This function returns a list that includes the response data, parameter estimates, and other data objects created by ‘ACER ConQuest’. These lists can be optionally coerced into R data frames.

# if no argument is provided to ConQuestSys, the example system file is read in by default.
myCqs <- ConQuestSys()
#> no system file provided, loading the example system file instead
#> System file is compressed. Decompressing ...
#> complete

You can see the data objects available within the object (e.g., str(myCqs)), and some useful objects will be:

gResponseData - The raw response data including information about the raw and key response to each item.
gYData - The regression data included in the estimation.
gAllCaseEstimates - the latent ability estimates (e.g., PVs, WLEs).
The option matrixout is available for several commands and will ensure data objects are saved to the system file - these can be found in gMatrixList. For example:
- The command matrixsampler and the option matrixout creates the objects matrixSampler_fit, matrixSampler_raw, matrix_userfit (the simulated data) in gMatrixList.
- The command itanal, and the option matrixout creates the objects *_counts, *_itemstats, *_ptbis, *_abilitymeansd (where “*” is a user defined prefix declared in the option matrixout) in gMatrixList .
- The command estimate and the option matrixout… for more information see the manual

Users can search the names of objects in the system file using the helper function searchConQuestSys:

# search for objects named history in myCqs
searchConQuestSys("history", myCqs)
#> [1] "gHistoryWeights" "gOldHistory"     "gBestHistory"    "gHistory"       
#> [5] "gCommandHistory"

Population statistics from plausible values

Large-scale assessments such as PISA report student achievement as a set of plausible values (PVs) and provide replicate weights for estimating sampling error. conquestr provides a family of functions that compute a statistic for each PV under the full-sample weight and under every replicate weight, then combine the results using Rubin’s rules:

pvStat() - the core function (any supported statistic, optionally by group)
pvMean(), pvSd(), pvVar(), pvMedian(), pvQuantile(), pvIqr() and pvDescribe()
pvDiff() - differences between groups (e.g., girls minus boys)
pvCor() - correlations between PVs and background variables
pvRegression() - weighted linear regression of PVs on background variables

Supported replication methods are "BRR", "Fay" (Fay’s variant of BRR, as used in PISA), "JK1" (delete-one jackknife) and "JK2" (paired jackknife, as used in TIMSS and PIRLS).

Simulate a PISA-style dataset

The example below simulates 64 schools in 32 variance strata (pairs of schools), 20 students per school, 10 PVs and Fay replicate weights (k = 0.5) built from a Hadamard matrix. PISA uses 80 replicates; 32 are used here to keep the example small.

set.seed(5309)
nStrata <- 32
nSchools <- nStrata * 2
nPerSchool <- 20
nPVs <- 10

# Sylvester construction of a 32 x 32 Hadamard matrix
hadamard <- matrix(1)
while (nrow(hadamard) < nStrata) hadamard <- rbind(cbind(hadamard, hadamard), cbind(hadamard, -hadamard))

school <- rep(seq_len(nSchools), each = nPerSchool)
stratum <- (school + 1) %/% 2
unit <- ifelse(school %% 2 == 1, 1, -1) # first or second school in the pair
n <- length(school)

pisaSim <- data.frame(
  school = school,
  sex = factor(sample(c("boy", "girl"), n, replace = TRUE)),
  escs = rnorm(nSchools, 0, 0.5)[school] + rnorm(n, 0, 0.8),
  W_FSTUWT = runif(n, 5, 50)
)
schoolEffect <- rnorm(nSchools, 0, 40)[school]
trueTheta <- 480 + 30 * pisaSim$escs + 12 * (pisaSim$sex == "girl") + schoolEffect + rnorm(n, 0, 70)
for (p in seq_len(nPVs)) {
  pisaSim[[paste0("PV", p, "READ")]] <- trueTheta + rnorm(n, 0, 25) # imputation noise
}

# Fay replicate weights: one school in each pair gets 1.5, the other 0.5
fayK <- 0.5
for (r in seq_len(nStrata)) {
  factorR <- ifelse(hadamard[stratum, r] * unit > 0, 2 - fayK, fayK)
  pisaSim[[paste0("W_FSTR", r)]] <- pisaSim$W_FSTUWT * factorR
}

pvs <- paste0("PV", seq_len(nPVs), "READ")
repWts <- paste0("W_FSTR", seq_len(nStrata))

Descriptive statistics

pvDescribe() returns the mean, SD, variance, quartiles and IQR. Each row reports the estimate, its standard error (se), the sampling variance U (samplingVar) and the between-PV variance B (imputationVar).

pvDescribe(pisaSim, pvs, weight = "W_FSTUWT", repWeights = repWts,
           repMethod = "Fay", fayK = fayK, by = "sex")
#>     sex  stat   estimate         se  samplingVar imputationVar nPVs nReps   n
#> 1   boy  mean  474.72213   4.615953 2.062039e+01  6.242144e-01   10    32 641
#> 2  girl  mean  492.81577   6.041746 3.521152e+01  1.173790e+00   10    32 639
#> 3   boy    sd   85.62484   3.084488 8.775039e+00  6.718418e-01   10    32 641
#> 4  girl    sd   91.17207   4.443683 1.899471e+01  6.832780e-01   10    32 639
#> 5   boy   var 7332.21824 524.164517 2.531133e+05  1.966830e+04   10    32 641
#> 6  girl   var 8312.96079 798.009731 6.119363e+05  2.262115e+04   10    32 639
#> 7   boy q0.25  417.07024   6.941046 3.946863e+01  7.917720e+00   10    32 641
#> 8  girl q0.25  432.98922   7.586118 5.056095e+01  6.352941e+00   10    32 639
#> 9   boy  q0.5  475.96979   5.755930 2.620271e+01  6.298201e+00   10    32 641
#> 10 girl  q0.5  493.03316   7.089010 4.478777e+01  4.969361e+00   10    32 639
#> 11  boy q0.75  531.27454   5.420923 2.443540e+01  4.500914e+00   10    32 641
#> 12 girl q0.75  549.38333   7.168611 4.739688e+01  3.629189e+00   10    32 639
#> 13  boy   iqr  114.20429   6.681075 3.344801e+01  1.017159e+01   10    32 641
#> 14 girl   iqr  116.39411   7.784710 5.017918e+01  9.475030e+00   10    32 639

imputationVar is B before Rubin’s correction, so the columns do not simply add up to se^2. The reported estimate is the average of only M PV estimates, and that average has its own simulation error with variance B / M. Rubin’s rules add it, so the total error variance is

\[\mathrm{SE}^2 = U + \left(1 + \tfrac{1}{M}\right) B.\]

With 10 PVs this adds 10% to B. To recover se from the columns:

desc <- pvDescribe(pisaSim, pvs, weight = "W_FSTUWT", repWeights = repWts,
                   repMethod = "Fay", fayK = fayK)
all.equal(desc$se, sqrt(desc$samplingVar + (1 + 1 / desc$nPVs) * desc$imputationVar))
#> [1] TRUE

The same holds for pvDiff(), pvCor() (on the z scale when combine = "fisherZ") and pvRegression(). See ?pvStat for details.

Group differences

pvDiff() computes the difference within each PV and replicate, so the standard error accounts for the covariance between groups that is induced by the sampling design.

pvDiff(pisaSim, pvs, by = "sex", pairs = list(c("girl", "boy")),
       weight = "W_FSTUWT", repWeights = repWts, repMethod = "Fay", fayK = fayK)
#>     contrast stat estimate       se samplingVar imputationVar        t
#> 1 girl - boy mean 18.09364 5.634755    29.94789      1.638702 3.211079

Correlations

pvCor() combines correlations on the raw r scale by default, or on the Fisher-z scale with combine = "fisherZ". The Fisher-z option gives confidence intervals that stay within [-1, 1].

pvCor(pisaSim, pvs, "escs", weight = "W_FSTUWT", repWeights = repWts,
      repMethod = "Fay", fayK = fayK, combine = "fisherZ")
#>   variable  estimate         se samplingVar imputationVar   ciLower   ciUpper
#> 1     escs 0.2801781 0.04637565 0.002369299   0.000148575 0.1870107 0.3683502
#>   zEstimate        zSe combine nPVs nReps    n
#> 1 0.2878754 0.05032625 fisherZ   10    32 1280

Regression

pvRegression() takes a one-sided formula of predictors; the PVs are the outcome.

pvRegression(pisaSim, pvs, ~ escs + sex, weight = "W_FSTUWT",
             repWeights = repWts, repMethod = "Fay", fayK = fayK)
#> 
#> Plausible-value regression
#> 
#> Call:
#> pvRegression(data = pisaSim, pvs = pvs, formula = ~escs + sex, 
#>     weight = "W_FSTUWT", repWeights = repWts, repMethod = "Fay", 
#>     fayK = fayK)
#> 
#> Coefficients:
#> (Intercept)         escs      sexgirl  
#>       476.8         27.1         16.8

Correlation and covariance matrices

pvCorMatrix() and pvCovMatrix() estimate the matrix for several measures at once. Each measure is either latent (its PV columns, in draw order) or manifest (a single observed column). Latent measures are paired by draw (PV1 with PV1, PV2 with PV2, …), which is appropriate when their PVs were drawn jointly from one multidimensional model. A manifest variable is the same in every draw, so its column is reused for each PV.

To illustrate, add a second latent domain (mathematics) drawn jointly with reading, and an observed variable (age in months):

trueMath <- 0.7 * (trueTheta - 480) + 470 + rnorm(n, 0, 60)
for (p in seq_len(nPVs)) {
  shared <- rnorm(n) # joint draw: reading and maths PV p share an error component
  pisaSim[[paste0("PV", p, "MATH")]] <- trueMath + 20 * shared + rnorm(n, 0, 15)
}
pisaSim$age <- round(rnorm(n, 189, 3.5))

measures <- list(
  read = paste0("PV", seq_len(nPVs), "READ"),
  math = paste0("PV", seq_len(nPVs), "MATH"),
  escs = "escs",
  age = "age"
)

The matrix format returns symmetric matrices of estimates, standard errors and confidence limits:

corMat <- pvCorMatrix(pisaSim, measures, weight = "W_FSTUWT", repWeights = repWts,
                      repMethod = "Fay", fayK = fayK, format = "matrix")
round(corMat$estimate, 3)
#>        read   math   escs    age
#> read  1.000  0.653  0.280 -0.029
#> math  0.653  1.000  0.165 -0.012
#> escs  0.280  0.165  1.000 -0.036
#> age  -0.029 -0.012 -0.036  1.000
round(corMat$se, 3)
#>       read  math  escs   age
#> read 0.000 0.027 0.047 0.029
#> math 0.027 0.000 0.049 0.029
#> escs 0.047 0.049 0.000 0.027
#> age  0.029 0.029 0.027 0.000

The long format has one row per pair. Pairs of two manifest variables are ordinary weighted correlations, so imputationVar is 0 and nPVs is 1:

pvCorMatrix(pisaSim, measures, weight = "W_FSTUWT", repWeights = repWts,
            repMethod = "Fay", fayK = fayK)[, c("dim1", "dim2", "estimate", "se",
                                                 "imputationVar", "nPVs")]
#>   dim1 dim2    estimate         se imputationVar nPVs
#> 1 read math  0.65324792 0.02732784  0.0001841724   10
#> 2 read escs  0.28014371 0.04678477  0.0001263433   10
#> 3 read  age -0.02853492 0.02922657  0.0000592743   10
#> 4 math escs  0.16526651 0.04878659  0.0000119018   10
#> 5 math  age -0.01200695 0.02889484  0.0001165340   10
#> 6 escs  age -0.03620490 0.02697831  0.0000000000    1

pvCovMatrix() returns variances on the diagonal; these match pvVar():

covMat <- pvCovMatrix(pisaSim, measures, weight = "W_FSTUWT", repWeights = repWts,
                      repMethod = "Fay", fayK = fayK, format = "matrix")
round(covMat$estimate, 1)
#>        read   math escs  age
#> read 7894.1 5151.4 22.7 -8.9
#> math 5151.4 7877.7 13.4 -3.7
#> escs   22.7   13.4  0.8 -0.1
#> age    -8.9   -3.7 -0.1 12.2
pvVar(pisaSim, measures$read, weight = "W_FSTUWT", repWeights = repWts,
      repMethod = "Fay", fayK = fayK)$estimate
#> [1] 7894.085

Two cautions apply when mixing latent and manifest measures:

A correlation between PVs and an observed variable is unbiased only if that variable (or a good proxy) was in the population (conditioning) model used to draw the PVs; otherwise it is biased towards zero. Here age was not used to generate the PVs and is unrelated to proficiency, so its correlations are near zero by construction. Check the study’s scaling documentation before interpreting such correlations.
With combine = "fisherZ", each element is back-transformed separately, so the matrix is not guaranteed to be positive semi-definite; a warning is given if it is not. The default, combine = "raw", always gives a valid matrix.