```
library(D3mirt)
library(mirt)
::knit_hooks$set(webgl = hook_webgl) knitr
```

`D3MIRT`

ModelingThe `D3mirt`

analysis is based on descriptive
multidimensional item response theory (DMIRT; Reckase, 2009, 1985;
Reckase & McKinley, 1991) and can be used to analyze dichotomous and
polytomous items (see Muraki & Carlson, 1995) in a three-dimensional
ability space. The method is foremost visual and illustrates item
characteristics with the help of vector geometry in which items are
represented by vector arrows.

In DMIRT analysis, also called within multidimensional modeling, it
is assumed that items in a multidimensional ability space can measure
single or multiple latent traits (Reckase, 2009, 1985; Reckase &
McKinley, 1991). The methodology is a type of data reduction technique
based on the compensatory model (Reckase, 2009), i.e., a type of
measurement model that uses linear combinations of \(\theta\)-values for ability assessment. The
method seeks to maximize item discrimination and so is
*descriptive* because the results describe the extent to which
items in a test are unidimensional, i.e., that the items discriminate on
one dimension only, or are within-multidimensional, i.e., that the items
discriminate on more than one dimension..

Regarding vector orientation, the angle of the vector arrows
indicates what traits, located along the orthogonal axes in the model,
an item can be said to describe (Reckase, 2009, 1985, Reckase &
McKinley, 1991). For instance, in a two-dimensional space, an item is
*unidimensional* if its item vector arrow is at \(0°\) with respect to one of the axes in the
model, and \(90°\) with respect to the
other. Such an item describes a singular trait only. In contrast, an
item is *within-multidimensional* if its item vector arrow is
oriented at \(45°\) in relation to the
axes in the model. Such an item describes both traits in the model
equally well. The same criteria are extended to the three-dimensional
case.

The DMIRT approach uses two types of item models, dependent on item
type. If dichotomous items are used, the analysis is based on the
two-parameter logistic model (M2PL). If polytomous items are used, the
analysis is extended to the two-parameter graded response model (MGRM;
Muraki & Carlson, 1995). In both cases, the estimation process
consists of first fitting a compensatory multidimensional two-parameter
model so that discrimination and difficulty (\(a\) and \(d\)) parameters can be extracted. For
`D3mirt`

, this implies that a compensatory three-dimensional
M2PL or MGRM must be fitted. Next, using these parameters the DMIRT
computation finds the direction of the highest possible discrimination
slope for each item. The direction, in turn, shows what latent trait the
item can be said to measure in the multidimensional latent model. The
output can be visualized as vector arrows representing item response
functions located in a multidimensional space pointing in the direction
of the maximum slope. The `D3mirt()`

function can be used to
calculate all necessary estimates for a three-dimensional DMIRT model,
as well as to plot the results.

The most central estimates, briefly explained below, in DMIRT analysis are the single multidimensional discrimination (\(MDISC\)) parameter, the multidimensional difficulty (\(MDIFF\)) parameters, and the directional discrimination (\(DDISC\)) parameters (Reckase2009, 1985; Reckase & McKinley, 1991). For reasons of simplicity, the fundamentals of DMIRT theory will be presented below limited to the two-dimensional case and the M2PL. The estimates will also be presented under the assumption that the axes \(l\) are orthogonal. This introduces the constraint that \(\sum_{l=1}^{m}cos^2\alpha_{jl} = 1\), i.e., that squared cosines equal \(1\) so that \(cos^2\, \alpha_{im} = 1-\sum_{k=1}^{m-1} cos^2\, \alpha_{ik}\), for item \(i\), person \(j\), and \(k\) orthogonal axes.

The \(MDISC\) for item \(i\) represents the highest level of discrimination the item \(i\) can achieve located in a multidimensional latent trait space (also called the point of inflection), with \(m\) number of dimensions and \(a_{ik}\) item slope parameters (Reckase, 2009, Reckase & McKinley, 1991).

\(MDISC = A_i = \sqrt {\sum_{k = 1}^{m} a^{2}_{ik} }\)

Similarly to unidimensional modeling, the \(MDISC\) parameter, also denoted as \(A_i\), is the parameter that set the steepness of the slope of the item response function at the point of inflection. The slope is, similarly to the unidimensional case, assessed multiplied with the constant \(\frac{1}{4}\) (omitted in the equation above).

The item orientation towards the steepest slope is set by taking the arc cosine of the ratio between the \(a_{il}\), i.e., the slope value of item \(i\) on the coordinate axis \(l\), and the \(MDISC\) (Reckase, 2009, 1985).

\[ a_{il}= cos^{-1}\left(\frac{a_{il}}{\sqrt{\sum_{k=1}^m a^2_{ik}}}\right) \]

The resulting direction vector is a characteristic of the item that describes the angular orientation of an item in a multidimensional latent trait space.

The multidimensional version of the difficulty parameter, here denoted \(B\) as the DMIRT counterpart to the \(b\)-parameter in the unidimensional item response theory model, for item \(i\), is defined as the negative intercept \(d_i\) divided by the \(MDISC\) (Reckase, 2009, 1985).

\[ MDIFF=B_i=\frac{-d_i}{\sqrt{\sum_{k=1}^m a^2_{ik}}} \]

The \(MDIFF\) is interpreted similarly as the difficulty parameter in the unidimensional model. That is, higher values indicate that higher levels of ability for a probability of a correct response of more than .5 are necessary. Moreover, the \(MDIFF\), just as in the unidimensional model, sets the distance from the origin of the model to the point of inflection. However, in DMIRT analysis, the \(MDIFF\) becomes a multidimensional location parameter that indicates the distance from the origin to the point of the steepest slope following the direction vector given by the direction vector equation.

For Likert items that hold multiple item response functions, the
\(MDIFF\) will be turned into an index
that indicates the distance from the origin to the point of inflection
for all response functions derived from an item. This, in turn, implies
that the \(MDIFF\) will show the
*difficulty range* for an item in the model.

The \(MDISC\) is visualized in DMIRT analysis by scaling the length of the vector arrows. In brief, the bottom location coordinates of the vector arrows, given by the \(MDIFF\) and direction vector, are multiplied with the \(MDISC\) so that items with higher \(MDISC\) have longer vector arrows (Reckase, 2009). Accordingly, shorter vector arrows indicate lower discrimination (which in turn indicates an increased amount of model violations of the item model).

Lastly, because the \(MDISC\) represents an item’s maximum level of discrimination in the model it is a global parameter that cannot be used to compare item discrimination locally. For the latter to be possible, the discrimination must be computed in a specified direction for the items. This can be achieved with the \(DDISC\), defined as follows (Reckase & McKinley, 1991).

\(DDISC = \sum_{k = 1}^{m} a_{ik} cos \alpha_{ik}\)

The \(DDISC\) describes the level of discrimination for one or more items in the model at the angle set by the direction of choice. Note, it is always true that \(DDISC \leq MDISC\).

In the `D3mirt`

package the \(DDISC\) is computed for all items following
the angle orientation given by the *construct vectors*.
Constructs, in this context, refer to the assumption that a subset of
items can measure a higher-order latent variable. In
`D3mirt`

, constructs are implemented as optional vectors
whose orientation is calculated as the average direction, i.e., the
average multidimensional discrimination of a subset of items in the
model. A construct vector will, therefore, point in the direction of the
maximum slope of an imaginary item response function indicated by the
items chosen by the user.

Item subsets for construct analysis are chosen by the user and can be exploratory (e.g., grouped based on observations) or theoretical (e.g., fixed a priori based on personality theory). If constructs are used, the output will include reporting of the \(DDISC\) parameter that shows how well the items discriminate under the assumption that they measure one of the constructs used in the analysis. That is, while the \(MDISC\) represents the maximum level of discrimination for the items in the model, the \(DDISC\) represents the local discrimination that makes it possible to compare item discrimination in a specific direction set by the constructs. Constructs are visually represented with construct vector arrows scaled to an arbitrary length.

The DMIRT method is currently limited to using the two-parameter
compensatory models. Regarding `D3mirt`

, the number of
dimensions is limited to three. However, the number of dimensions need
not be exactly three, but up to three. This since only two items in the
set are needed to identify the model (see the section on model
identification below for more on the statistical requirements).

The package includes the following functions.

`modid()`

: D3mirt Model Identification`D3mirt()`

: 3D DMIRT Model Estimation`summary()`

: Summary Function for`D3mirt()`

`plotD3mirt()`

: Graphical Output for`D3mirt()`

In what follows, the `D3mirt`

procedure and some of its
functions and options will be described using the built-in data set
“anes0809offwaves”. The data set (\(N = 1046,
M_{age} = 51.33, SD = 14.56, 57\%\) Female) is a subset from the
American National Election Survey (ANES) from the 2008-2009 Panel Study
Off Wave Questionnaires, December 2009 (DeBell, et al, 2010; https://electionstudies.org/data-center/2008-2009-panel-study/).
All items measure moral preferences and are positively scored of Likert
type, ranging from 1 = *Strongly Disagree* to 6 = *Strongly
Agree*. Demographic variables include age and gender
(male/female).

The sections below are sorted under the following headings.

- Model Identification

- 1.1 The
`modid()`

Function- Step 1: Explore the Data Structure
- Step 2: Item Selection

- 1.2 The Model Identification Procedure
- 1.3 Criteria For Model Identification
- 1.4 Limitations
- 1.3 Item Selection

`D3mirt`

Model Estimation

- 2.1 The
`D3mirt()`

Function- Fitting the Compensatory Multidimensional Model
- The
`D3mirt`

Function Call - Constructs

- Plotting

- 3.1 The
`plotD3mirt()`

Function- Unidimensionality vs. Within-Dimensionality
- Model Violations
- Illustration: A Short Item and Dimensionality Analysis
- Other Graphical Options
`items`

`scale`

`D3mirt`

Profile Analysis

`modid()`

FunctionFor the `D3mirt`

analysis, two items with the following
properties must be chosen to identify the compensatory model (Reckase,
2009).

The first item should not load on the second and third axes (\(y\) and \(z\)), while the second item should not load
on the third axis (\(z\)). If this can
be empirically achieved, it will be possible to create the orthogonal
structure necessary for the analysis. The `modid()`

function
can help by suggesting what items to use for this purpose.

Note, if improper items are chosen the model will be hard to
interpret, or even unempirical if the data is forced on the model by
compromising the necessary statistical requirements. The
`modid()`

function was, therefore, developed to hinder the
latter and to maximize the interpretive meaning by using an algorithmic
approach. Briefly, the function first order factors by the sum of
squares and from this select the strongest loading items that meet the
statistical assumptions described above. This orders the model so that
the strongest loading item, from the strongest factor, always aligns
perfectly with the \(x\)-axis and that
the other items follow thereon. The model identification process is
described in more detail below.

If the model is not known, the factor structure must be explored with
exploratory factor analysis (EFA). However, because `D3mirt`

analysis is based on the M2PL and the MGRM, it is recommended to use
multidimensional item response theory EFA methods, such as the EFA
option in `mirt::mirt`

(Chalmers, 2012) with
`ìtemtype = 'graded'`

or `'2PL'`

, so that the EFA
is performed with the proper item model. Note, the EFA is only used to
find model identification items that meet the necessary DMIRT model
specification requirements. The EFA model itself is discarded after this
step in the procedure. It user can therefore try different types of
rotation methods and compare model identification results. Note, because
the EFA in the `mirt()`

function takes a long time to
perform, the item loadings from the EFA for this example are stored in
the package file “efa.Rdata” and can be loaded in the next code chunk
below.

```
# Load data
data("anes0809offwaves")
<- anes0809offwaves
x <- x[,3:22] # Remove columns for age and gender
x
# Fit a three-factor EFA model with the mirt package
<- mirt::mirt(x, 3, itemtype = 'graded')
e
# Assign data frame with factor loadings with oblimin rotation
<- summary(e, rotate= 'oblimin')
f <- data.frame(f$rotF) h
```

The `modid()`

takes in the factor solution from the EFA,
assigned to a data frame \(h\), and
outputs an \(S3\) object of class
`modid`

containing lists with data frames of estimates. The
most important is the item lists, (denoted `$id`

), containing
data frames that present suggestions on what items
(`item.1`

…`item.n`

) to use for the model
identification. The data frames have one column for the loadings from
each item on the factor of interest, and one column with absolute sum
scores (denoted `ABS`

) for each item calculated from the
remaining factor loadings in the model. Each item list is sorted with
the lowest absolute sum score highest up. Consequently, the top items in
each list are the items that best meet the necessary statistical
requirements for model identification. Therefore, for a
three-dimensional model, all else equal, the item highest up in the
first list should be used to identify the \(x\)-axis, and the item highest up in the
second list should be used to identify the \(y\)-axis, and so on.

```
# Optional: Load the EFA data for this example directly from the package file
load("efa.Rdata")
# Call to modid()
modid(h)
#> $id
#> $id[[1]]
#> Item.1 ABS
#> W7Q3 0.8490903 0.01267429
#> W7Q5 0.8068828 0.04904288
#> W7Q1 0.7543773 0.07203403
#> W7Q2 0.8727652 0.09042323
#>
#> $id[[2]]
#> Item.2 ABS
#> W7Q20 0.7858844 0.000653665
#> W7Q18 0.6812006 0.080729410
#>
#>
#> $ss.loadings
#> F1 F3 F2
#> 5.328901 2.113441 1.689176
#>
#> $loadings
#> F1 F3 F2
#> W7Q11 0.21764339 0.101650054 0.535617304
#> W7Q12 0.07445495 -0.079798063 0.554858770
#> W7Q13 -0.01369236 -0.018393288 0.769725868
#> W7Q14 -0.03993821 0.145349221 0.564359537
#> W7Q15 0.10245563 0.453634967 -0.099027661
#> W7Q16 0.16609422 0.212788638 0.126237569
#> W7Q17 0.21251128 0.576133340 0.039833393
#> W7Q18 -0.05188854 0.681200616 -0.080729410
#> W7Q19 0.02592854 0.626381734 0.125087323
#> W7Q20 -0.05079509 0.785884397 0.000653665
#> W7Q1 0.75437734 0.040367304 -0.031666723
#> W7Q2 0.87276522 -0.024432875 -0.065990357
#> W7Q3 0.84909025 -0.010993866 0.001680422
#> W7Q4 0.66228706 0.032788311 0.101713685
#> W7Q5 0.80688278 -0.040279704 -0.008763174
#> W7Q6 0.66856685 0.054813498 0.102271288
#> W7Q7 0.56078396 -0.013762611 0.211076266
#> W7Q8 0.56779935 0.042979814 0.204500105
#> W7Q9 0.60483387 0.090013632 0.088259630
#> W7Q10 0.77064478 0.009554713 -0.116375618
#>
#> attr(,"class")
#> [1] "modid"
```

As can be seen, the first item, that will be used to identify the
\(x\)-axis, is found in the first list,
`id[[1]]`

and `item.1`

. In this case, the best
item for the \(x\)-axis is item “W7Q3”.
The item that identifies the \(y\)-axis
is found in the next list, `id[[2]]`

and `Item.2`

.
In this case, the best item for the \(y\)-axis is item “W7Q20”.

Sometimes, however, the model is hard to identify. If this happens, try the following in the order suggested. For more on the function arguments see the section below regarding the model identification procedure.

- Change the rotation method in the EFA, e.g., to change from
*oblimin*to*varimax*. - Adjust the
`lower`

bound in`modid()`

. - Override factor order with
`fac.order`

. - Adjust the
`upper`

bound in`modid()`

.

The latter (point 4) should, however, only be used as a last resort. This is because the upper bound sets the upper limit for item inclusion. Adjusting this limit too high means that the necessary statistical requirements are compromised. The lower limit (point 2), however, only increases the size of the item pool used for the item selection. It is, therefore, recommended to adjust the lower limit up and down to see if the output differs, and from that make the final decision.

```
# Call to modid with increased lower and higher bound
modid(h, lower = 1, upper = .12 )
#> $id
#> $id[[1]]
#> Item.1 ABS
#> W7Q3 0.8490903 0.01267429
#> W7Q5 0.8068828 0.04904288
#> W7Q1 0.7543773 0.07203403
#> W7Q2 0.8727652 0.09042323
#>
#> $id[[2]]
#> Item.2 ABS
#> W7Q20 0.7858844 0.000653665
#> W7Q17 0.5761333 0.039833393
#> W7Q18 0.6812006 0.080729410
#> W7Q15 0.4536350 0.099027661
#>
#>
#> $ss.loadings
#> F1 F3 F2
#> 5.328901 2.113441 1.689176
#>
#> $loadings
#> F1 F3 F2
#> W7Q11 0.21764339 0.101650054 0.535617304
#> W7Q12 0.07445495 -0.079798063 0.554858770
#> W7Q13 -0.01369236 -0.018393288 0.769725868
#> W7Q14 -0.03993821 0.145349221 0.564359537
#> W7Q15 0.10245563 0.453634967 -0.099027661
#> W7Q16 0.16609422 0.212788638 0.126237569
#> W7Q17 0.21251128 0.576133340 0.039833393
#> W7Q18 -0.05188854 0.681200616 -0.080729410
#> W7Q19 0.02592854 0.626381734 0.125087323
#> W7Q20 -0.05079509 0.785884397 0.000653665
#> W7Q1 0.75437734 0.040367304 -0.031666723
#> W7Q2 0.87276522 -0.024432875 -0.065990357
#> W7Q3 0.84909025 -0.010993866 0.001680422
#> W7Q4 0.66228706 0.032788311 0.101713685
#> W7Q5 0.80688278 -0.040279704 -0.008763174
#> W7Q6 0.66856685 0.054813498 0.102271288
#> W7Q7 0.56078396 -0.013762611 0.211076266
#> W7Q8 0.56779935 0.042979814 0.204500105
#> W7Q9 0.60483387 0.090013632 0.088259630
#> W7Q10 0.77064478 0.009554713 -0.116375618
#>
#> attr(,"class")
#> [1] "modid"
```

In this case, the same items are suggested even after adjusting both
the `lower`

and `upper`

bounds.

Another option (point 3) is to override the factor order with the
`fac.order`

argument. More specifically, `modid`

orders factor by squared factor loadings so that the strongest factor is
used first, the second strongest factor is used second, and so on.
Sometimes, however, there is only a very small difference between the
squared factor loadings that in turn can cause problems (often only
observable at later stages) when trying to find the best items for the
model identification. In such situations, it can be useful to rearrange
the factor order manually to see if the model solution improves.

Since the squared factor loadings for factors 2 and 3 in this example are somewhat similar, it could be useful to compare the results. The latter is, however, outside of the scope of this vignette.

```
# Override factor order by reversing columns in the original data frame
modid(h, fac.order = c(3, 2, 1))
#> $id
#> $id[[1]]
#> Item.1 ABS
#> W7Q20 0.7858844 0.05144876
#>
#> $id[[2]]
#> Item.2 ABS
#> W7Q13 0.7697259 0.01369236
#> W7Q14 0.5643595 0.03993821
#>
#>
#> $ss.loadings
#> F3 F2 F1
#> 2.113441 1.689176 5.328901
#>
#> $loadings
#> F3 F2 F1
#> W7Q11 0.101650054 0.535617304 0.21764339
#> W7Q12 -0.079798063 0.554858770 0.07445495
#> W7Q13 -0.018393288 0.769725868 -0.01369236
#> W7Q14 0.145349221 0.564359537 -0.03993821
#> W7Q15 0.453634967 -0.099027661 0.10245563
#> W7Q16 0.212788638 0.126237569 0.16609422
#> W7Q17 0.576133340 0.039833393 0.21251128
#> W7Q18 0.681200616 -0.080729410 -0.05188854
#> W7Q19 0.626381734 0.125087323 0.02592854
#> W7Q20 0.785884397 0.000653665 -0.05079509
#> W7Q1 0.040367304 -0.031666723 0.75437734
#> W7Q2 -0.024432875 -0.065990357 0.87276522
#> W7Q3 -0.010993866 0.001680422 0.84909025
#> W7Q4 0.032788311 0.101713685 0.66228706
#> W7Q5 -0.040279704 -0.008763174 0.80688278
#> W7Q6 0.054813498 0.102271288 0.66856685
#> W7Q7 -0.013762611 0.211076266 0.56078396
#> W7Q8 0.042979814 0.204500105 0.56779935
#> W7Q9 0.090013632 0.088259630 0.60483387
#> W7Q10 0.009554713 -0.116375618 0.77064478
#>
#> attr(,"class")
#> [1] "modid"
```

In this case, we find the item that was previously suggested for the \(y\)-axis, is now suggested for the \(x\)-axis. For the \(y\)-axis, we find a new item, “W7Q13”, suggested to identify the \(y\)-axis.

The `modid()`

function uses an iterating algorithmic
procedure that can be user adjusted. In brief, in the default automatic
mode, `modid()`

starts by first calculating the sum of
squares loadings on all factors \(F\)
in the data frame \(x\) and then
rearranging the columns in \(x\), in
decreasing order following the level of strength of the sum of squares
loadings. Next, the function creates a list containing the factor
loadings on the first factor, \(f_1\),
and absolute sum scores of the factor loadings on the remaining factors,
i.e., \(F-f_1\), row-wise. The list is
rearranged in decreasing order based row-wise on factor loading strength
on \(f_1\).

Next, items are selected by first scaling \(f_1\) loadings and then extracting the items with the highest loadings on \(f_1\), up to a standard deviation of \(0.5\) (the default setting) as the lower bound criteria counting from the top. That is, rows with raw factor scores and absolute sum scores are extracted until the lower bound for factor loadings on \(f_1\) is reached. This allows the function to extract more rows in the case empirical factor loadings are very similar in strength. The result is recorded as a nested list before the function starts over with the next factor, \(f_2\), and so on.

For every iteration, the algorithm jumps to the next factor in the EFA model, rearranges rows, and extracts the strongest loading items. However the absolute sum score is always assessed based on the number of factors less than the total number of factors, following the order of iteration, That is, iteration \(1\) use factor loadings from all factors \(F-f_1\), iteration \(2\) \(F-(f_{1,2})\), iteration \(3\) \(F-(f_{1,2,3})\), and so on, when calculating the absolute sum scores.

Model identification items should preferably (a) have an absolute sum score of \(\leq .10\) and (b) have the highest factor loading scores on the factor of interest. Of these two criteria, (a) should be given the strongest weight in the selection decision. If these conditions cannot be met, the user is advised to proceed with caution since the loading scores, therefore, imply that an adequate orthogonal structure may not be empirically attainable. If problems in the model identification process occur, please follow the advice given above.

The `modid()`

function is not limited to three-dimensional
analysis and can be used on any number of dimensions. Although based on
suggestions on model identification given by Reckase (2009) for this
type of analysis, the function offers some expansions that introduce
more precision. The latter foremost consists in incorporating the sum of
squares and factor loadings in the item selection process (unless the
user has not specified otherwise). Experience tells us that this is good
practice that often leads to better results compared to other known
options. However, it is important to recognize that the model
identification procedure only gives suggestions to the model
specification, and there could be situations where the researcher should
consider other methods. Note, that two items can be found to identify
the model do not imply successful outcomes when using this methodology
(i.e., that the model is *good*). But it does suggest that the
methodology *can* be used and the results will be possible to
interpret in a meaningful way.

`D3mirt()`

FunctionThe `D3mirt()`

function takes in a data frame with model
parameters from a three-dimensional compensatory model and returns an
\(S3\) object of class
`D3mirt`

with lists of \(a\)
and \(d\), \(MDISC\), and \(MDIFF\) parameters, direction cosines, and
spherical coordinates. Regarding the latter, spherical coordinates are
represented by \(\theta\) and \(\phi\). The \(\theta\) coordinate is the positive or
negative angle in degrees, starting from the \(x\)-axis, of the vector projections from
the vector arrows in the \(xz\)-plane
up to \(\pm 180°\). Note, the \(\theta\) angle is oriented following the
positive pole of the \(x\) and \(z\) axis so that the angle increases
clockwise in the graphical output. The \(\phi\) coordinate is the positive angle in
degrees from the \(y\)-axis and the
vectors. Note, the \(\rho\) coordinate
from the spherical coordinate system is in DMIRT represented by the
\(MDIFF\), and so is reported
separately.

If constructs are used, the function also returns construct direction cosines, spherical coordinates for the construct vector arrows, and \(DDISC\) parameters (one index per construct).

The three-dimensional compensatory model is specified so that all
items load on all three factors in the model, and all factors are
constrained to be orthogonal (see below). The fitting of the model is
preferably done with the `mirt::mirt`

(Chalmers, 2012)
function. Observe that the START and FIXED commands are used to fix the
slope parameters on the second, \(a2\),
and third, \(a3\). factor for item
\(_1\) (W7Q3), and the slope on the
third, \(a3\), factor for item \(i_2\) (W7Q20). Note, because the fitting of
the compensatory model in the `mirt()`

function takes a long
time to fit, the item parameters are available in the package file
“d.Rdata” and can be loaded in the next code chunk below.

```
# Load data
data("anes0809offwaves")
<- anes0809offwaves
x <- x[,3:22] # Remove columns for age and gender
x
# Fit a three-dimensional graded response model with orthogonal factors
# Example below uses Likert items from the built-in data set "anes0809offwaves"
# Item W7Q3 and item W7Q20 was selected with `modid()`
# The model specification set all items in the data set (1-20)
# to load on all three factors (F1-F3)
# The START and FIXED commands are used on the two items to identify the DMIRT model
<- ' F1 = 1-20
spec F2 = 1-20
F3 = 1-20
START=(W7Q3,a2,0)
START=(W7Q3,a3,0)
START=(W7Q20,a3,0)
FIXED=(W7Q3,a2)
FIXED=(W7Q3,a3)
FIXED=(W7Q20,a3) '
<- mirt::mirt(x,
mod1
spec, itemtype = 'graded',
SE = TRUE,
method = 'QMCEM')
# Assign a data frame with factor loadings (located in the first three columns in mod1),
# and difficulty parameters (columns 4-8 in mod1) with mirt::coef and $'items'[,1:8]))
<- data.frame(mirt::coef(mod1,
d simplify=TRUE)$'items'[,1:8])
```

`D3mirt()`

Function CallThe `D3mirt()`

function call is straightforward. The
output, however, is lengthy so the use of the generic summary function,
included in the package, when inspecting the results is recommended.

```
# Optional: Load data frame d for for this example directly from the package file
load("d.Rdata")
# Call D3mirt()
<- D3mirt(d)
g summary(g) # Show summary of results
#>
#> 3Dmirt object with 20 items and 5 levels of difficulty
#> $model.est
#> a1 a2 a3 d1 d2 d3 d4 d5
#> W7Q11 1.4237 0.4675 1.0439 6.2189 4.6926 3.5435 1.1920 -1.8576
#> W7Q12 0.7604 0.0410 0.9367 4.1360 2.8771 2.3419 1.1790 -0.4240
#> W7Q13 1.1278 0.2911 1.6930 5.8892 4.3988 3.4413 1.8946 -0.6008
#> W7Q14 0.7447 0.4829 0.9785 5.3891 3.9333 3.0258 0.8143 -1.5868
#> W7Q15 0.4551 0.7870 -0.1606 4.3207 3.0545 2.3969 0.9187 -0.9705
#> W7Q16 0.6237 0.4140 0.1798 3.7249 2.0305 1.1658 -0.0612 -1.8085
#> W7Q17 1.1892 1.3412 0.0563 6.9013 5.8023 4.9345 2.7916 -0.0041
#> W7Q18 0.4106 1.3542 -0.1369 3.7837 2.0985 1.4183 0.1828 -1.9855
#> W7Q19 0.8580 1.4099 0.2279 4.4978 2.6483 1.6730 0.3741 -1.9966
#> W7Q20 0.7357 1.9067 0.0000 4.6378 2.3633 1.2791 -0.3431 -2.9190
#> W7Q1 2.0298 0.1643 -0.1231 8.0865 7.0640 5.9876 3.2015 -0.4835
#> W7Q2 2.6215 -0.0027 -0.2582 9.2885 6.6187 4.5102 1.6649 -2.4439
#> W7Q3 2.7923 0.0000 0.0000 10.4894 7.5887 5.6776 2.7172 -1.1789
#> W7Q4 1.9045 0.1875 0.1495 7.3750 6.0465 4.9813 2.4830 -1.1145
#> W7Q5 2.2425 -0.0287 -0.0839 8.4279 6.6712 4.9049 1.8253 -1.8316
#> W7Q6 2.0021 0.2390 0.1571 8.0684 6.3577 4.9520 2.3300 -1.0189
#> W7Q7 1.6286 0.1034 0.3595 6.0178 4.8974 3.6908 1.6326 -1.3484
#> W7Q8 1.7774 0.2252 0.3531 6.9171 5.1822 3.7661 1.4844 -1.8332
#> W7Q9 1.7198 0.2494 0.1281 7.5586 4.9755 3.3648 0.9343 -2.2094
#> W7Q10 1.7696 0.1272 -0.1406 8.3638 5.7397 4.2863 1.9647 -0.6642
#>
#> $dmirt.est
#> MDISC MDIFF1 MDIFF2 MDIFF3 MDIFF4 MDIFF5
#> W7Q11 1.8263 -3.4052 -2.5695 -1.9403 -0.6527 1.0171
#> W7Q12 1.2072 -3.4262 -2.3833 -1.9400 -0.9767 0.3512
#> W7Q13 2.0550 -2.8658 -2.1405 -1.6746 -0.9220 0.2923
#> W7Q14 1.3211 -4.0793 -2.9773 -2.2904 -0.6164 1.2011
#> W7Q15 0.9232 -4.6801 -3.3085 -2.5963 -0.9951 1.0512
#> W7Q16 0.7699 -4.8382 -2.6374 -1.5142 0.0795 2.3490
#> W7Q17 1.7934 -3.8481 -3.2353 -2.7514 -1.5566 0.0023
#> W7Q18 1.4217 -2.6613 -1.4760 -0.9976 -0.1286 1.3966
#> W7Q19 1.6661 -2.6996 -1.5896 -1.0042 -0.2245 1.1984
#> W7Q20 2.0437 -2.2693 -1.1564 -0.6259 0.1679 1.4283
#> W7Q1 2.0402 -3.9637 -3.4625 -2.9348 -1.5692 0.2370
#> W7Q2 2.6342 -3.5261 -2.5126 -1.7122 -0.6320 0.9278
#> W7Q3 2.7923 -3.7565 -2.7177 -2.0333 -0.9731 0.4222
#> W7Q4 1.9195 -3.8421 -3.1500 -2.5950 -1.2935 0.5806
#> W7Q5 2.2442 -3.7554 -2.9726 -2.1856 -0.8133 0.8161
#> W7Q6 2.0225 -3.9894 -3.1436 -2.4485 -1.1521 0.5038
#> W7Q7 1.6710 -3.6013 -2.9308 -2.2087 -0.9770 0.8069
#> W7Q8 1.8261 -3.7880 -2.8379 -2.0624 -0.8129 1.0039
#> W7Q9 1.7425 -4.3377 -2.8553 -1.9310 -0.5362 1.2679
#> W7Q10 1.7798 -4.6994 -3.2250 -2.4083 -1.1039 0.3732
#>
#> $dmirt.angles
#> D.Cos X D.Cos Y D.Cos Z Theta Phi
#> W7Q11 0.7796 0.2560 0.5716 36.2490 75.1671
#> W7Q12 0.6299 0.0340 0.7759 50.9315 88.0517
#> W7Q13 0.5488 0.1417 0.8238 56.3293 81.8559
#> W7Q14 0.5637 0.3656 0.7407 52.7295 68.5585
#> W7Q15 0.4929 0.8525 -0.1739 -19.4363 31.5157
#> W7Q16 0.8101 0.5377 0.2336 16.0827 57.4726
#> W7Q17 0.6631 0.7479 0.0314 2.7110 41.5948
#> W7Q18 0.2888 0.9525 -0.0963 -18.4366 17.7258
#> W7Q19 0.5150 0.8462 0.1368 14.8745 32.1959
#> W7Q20 0.3600 0.9330 0.0000 0.0000 21.0985
#> W7Q1 0.9949 0.0806 -0.0603 -3.4700 85.3796
#> W7Q2 0.9952 -0.0010 -0.0980 -5.6253 90.0588
#> W7Q3 1.0000 0.0000 0.0000 0.0000 90.0000
#> W7Q4 0.9922 0.0977 0.0779 4.4875 84.3945
#> W7Q5 0.9992 -0.0128 -0.0374 -2.1436 90.7317
#> W7Q6 0.9899 0.1182 0.0777 4.4863 83.2121
#> W7Q7 0.9746 0.0619 0.2152 12.4486 86.4516
#> W7Q8 0.9733 0.1233 0.1933 11.2350 82.9150
#> W7Q9 0.9870 0.1431 0.0735 4.2594 81.7720
#> W7Q10 0.9943 0.0715 -0.0790 -4.5437 85.9000
```

As can be seen, the output from the summary function starts by reporting the number of items and difficulty levels. If constructs were used, the output will also include statements on what items belong to which construct (see example below). Next, the factor loadings and the difficulty parameters from the compensatory model are reported followed by all necessary DMIRT estimates. Examples of how DMIRT estimates can be used when reporting results are given below in the item and dimensionality analysis example.

Constructs can be included in the analysis by creating one or more
nested lists that indicate what items belong to what construct. From
this, the `D3mirt()`

function finds the average direction of
the subset of items contained in each construct list, by adding and
normalizing the direction cosines for the items and scaling the
construct direction vector to an arbitrary length. Note, the length of
the construct vector arrows can be user adjusted.

The construct vector arrows can contribute to the analysis by (a) visualizing the average direction for a subset set of items, and (b) showing how all items discriminate locally in the direction of the construct vector with the help of the \(DDISC\) index.

The constructs included below were grouped based on exploratory reasons, i.e., because these items cluster in the model (will be observable in the graphical output).

```
# Call to D3mirt(), including optional nested lists for three constructs
# Item W7Q16 is not included in any construct because of model violations
# The model violations for the W7Q16 item can be seen when plotting the model
<- list(list (1,2,3,4),
c list(5,7,8,9,10),
list(11,12,13,14,15,15,16,17,18,19,20))
<- D3mirt(d, c)
g summary(g)
#>
#> 3Dmirt object with 20 items and 5 levels of difficulty
#>
#> Construct vector 1 contains items 1, 2, 3, 4
#>
#> Construct vector 2 contains items 5, 7, 8, 9, 10
#>
#> Construct vector 3 contains items 11, 12, 13, 14, 15, 15, 16, 17, 18, 19, 20
#> $model.est
#> a1 a2 a3 d1 d2 d3 d4 d5
#> W7Q11 1.4237 0.4675 1.0439 6.2189 4.6926 3.5435 1.1920 -1.8576
#> W7Q12 0.7604 0.0410 0.9367 4.1360 2.8771 2.3419 1.1790 -0.4240
#> W7Q13 1.1278 0.2911 1.6930 5.8892 4.3988 3.4413 1.8946 -0.6008
#> W7Q14 0.7447 0.4829 0.9785 5.3891 3.9333 3.0258 0.8143 -1.5868
#> W7Q15 0.4551 0.7870 -0.1606 4.3207 3.0545 2.3969 0.9187 -0.9705
#> W7Q16 0.6237 0.4140 0.1798 3.7249 2.0305 1.1658 -0.0612 -1.8085
#> W7Q17 1.1892 1.3412 0.0563 6.9013 5.8023 4.9345 2.7916 -0.0041
#> W7Q18 0.4106 1.3542 -0.1369 3.7837 2.0985 1.4183 0.1828 -1.9855
#> W7Q19 0.8580 1.4099 0.2279 4.4978 2.6483 1.6730 0.3741 -1.9966
#> W7Q20 0.7357 1.9067 0.0000 4.6378 2.3633 1.2791 -0.3431 -2.9190
#> W7Q1 2.0298 0.1643 -0.1231 8.0865 7.0640 5.9876 3.2015 -0.4835
#> W7Q2 2.6215 -0.0027 -0.2582 9.2885 6.6187 4.5102 1.6649 -2.4439
#> W7Q3 2.7923 0.0000 0.0000 10.4894 7.5887 5.6776 2.7172 -1.1789
#> W7Q4 1.9045 0.1875 0.1495 7.3750 6.0465 4.9813 2.4830 -1.1145
#> W7Q5 2.2425 -0.0287 -0.0839 8.4279 6.6712 4.9049 1.8253 -1.8316
#> W7Q6 2.0021 0.2390 0.1571 8.0684 6.3577 4.9520 2.3300 -1.0189
#> W7Q7 1.6286 0.1034 0.3595 6.0178 4.8974 3.6908 1.6326 -1.3484
#> W7Q8 1.7774 0.2252 0.3531 6.9171 5.1822 3.7661 1.4844 -1.8332
#> W7Q9 1.7198 0.2494 0.1281 7.5586 4.9755 3.3648 0.9343 -2.2094
#> W7Q10 1.7696 0.1272 -0.1406 8.3638 5.7397 4.2863 1.9647 -0.6642
#>
#> $dmirt.est
#> MDISC MDIFF1 MDIFF2 MDIFF3 MDIFF4 MDIFF5
#> W7Q11 1.8263 -3.4052 -2.5695 -1.9403 -0.6527 1.0171
#> W7Q12 1.2072 -3.4262 -2.3833 -1.9400 -0.9767 0.3512
#> W7Q13 2.0550 -2.8658 -2.1405 -1.6746 -0.9220 0.2923
#> W7Q14 1.3211 -4.0793 -2.9773 -2.2904 -0.6164 1.2011
#> W7Q15 0.9232 -4.6801 -3.3085 -2.5963 -0.9951 1.0512
#> W7Q16 0.7699 -4.8382 -2.6374 -1.5142 0.0795 2.3490
#> W7Q17 1.7934 -3.8481 -3.2353 -2.7514 -1.5566 0.0023
#> W7Q18 1.4217 -2.6613 -1.4760 -0.9976 -0.1286 1.3966
#> W7Q19 1.6661 -2.6996 -1.5896 -1.0042 -0.2245 1.1984
#> W7Q20 2.0437 -2.2693 -1.1564 -0.6259 0.1679 1.4283
#> W7Q1 2.0402 -3.9637 -3.4625 -2.9348 -1.5692 0.2370
#> W7Q2 2.6342 -3.5261 -2.5126 -1.7122 -0.6320 0.9278
#> W7Q3 2.7923 -3.7565 -2.7177 -2.0333 -0.9731 0.4222
#> W7Q4 1.9195 -3.8421 -3.1500 -2.5950 -1.2935 0.5806
#> W7Q5 2.2442 -3.7554 -2.9726 -2.1856 -0.8133 0.8161
#> W7Q6 2.0225 -3.9894 -3.1436 -2.4485 -1.1521 0.5038
#> W7Q7 1.6710 -3.6013 -2.9308 -2.2087 -0.9770 0.8069
#> W7Q8 1.8261 -3.7880 -2.8379 -2.0624 -0.8129 1.0039
#> W7Q9 1.7425 -4.3377 -2.8553 -1.9310 -0.5362 1.2679
#> W7Q10 1.7798 -4.6994 -3.2250 -2.4083 -1.1039 0.3732
#>
#> $dmirt.angles
#> D.Cos X D.Cos Y D.Cos Z Theta Phi
#> W7Q11 0.7796 0.2560 0.5716 36.2490 75.1671
#> W7Q12 0.6299 0.0340 0.7759 50.9315 88.0517
#> W7Q13 0.5488 0.1417 0.8238 56.3293 81.8559
#> W7Q14 0.5637 0.3656 0.7407 52.7295 68.5585
#> W7Q15 0.4929 0.8525 -0.1739 -19.4363 31.5157
#> W7Q16 0.8101 0.5377 0.2336 16.0827 57.4726
#> W7Q17 0.6631 0.7479 0.0314 2.7110 41.5948
#> W7Q18 0.2888 0.9525 -0.0963 -18.4366 17.7258
#> W7Q19 0.5150 0.8462 0.1368 14.8745 32.1959
#> W7Q20 0.3600 0.9330 0.0000 0.0000 21.0985
#> W7Q1 0.9949 0.0806 -0.0603 -3.4700 85.3796
#> W7Q2 0.9952 -0.0010 -0.0980 -5.6253 90.0588
#> W7Q3 1.0000 0.0000 0.0000 0.0000 90.0000
#> W7Q4 0.9922 0.0977 0.0779 4.4875 84.3945
#> W7Q5 0.9992 -0.0128 -0.0374 -2.1436 90.7317
#> W7Q6 0.9899 0.1182 0.0777 4.4863 83.2121
#> W7Q7 0.9746 0.0619 0.2152 12.4486 86.4516
#> W7Q8 0.9733 0.1233 0.1933 11.2350 82.9150
#> W7Q9 0.9870 0.1431 0.0735 4.2594 81.7720
#> W7Q10 0.9943 0.0715 -0.0790 -4.5437 85.9000
#>
#> $construct.angles
#> C.Cos X C.Cos Y C.Cos Z Theta Phi
#> C1 0.6411 0.2026 0.7402 49.1064 78.3081
#> C2 0.4720 0.8814 -0.0208 -2.5190 28.1921
#> C3 0.9977 0.0613 0.0298 1.7098 86.4857
#>
#> $ddisc
#> DDISC1 DDISC2 DDISC3
#> W7Q11 1.7802 1.0624 1.4802
#> W7Q12 1.1892 0.3756 0.7890
#> W7Q13 2.0353 0.7537 1.1935
#> W7Q14 1.2996 0.7568 0.8017
#> W7Q15 0.3324 0.9118 0.4975
#> W7Q16 0.6169 0.6555 0.6530
#> W7Q17 1.0759 1.7422 1.2704
#> W7Q18 0.4364 1.3902 0.4886
#> W7Q19 1.0044 1.6428 0.9492
#> W7Q20 0.8580 2.0277 0.8509
#> W7Q1 1.2434 1.1054 2.0315
#> W7Q2 1.4889 1.2403 2.6076
#> W7Q3 1.7901 1.3179 2.7858
#> W7Q4 1.3696 1.0610 1.9160
#> W7Q5 1.3697 1.0349 2.2330
#> W7Q6 1.4482 1.1524 2.0168
#> W7Q7 1.3311 0.8523 1.6418
#> W7Q8 1.4464 1.0301 1.7976
#> W7Q9 1.2479 1.0288 1.7349
#> W7Q10 1.0561 0.9503 1.7691
```

Compared to the previous function call to `D3mirt()`

, we
now see two extra frames: one frame with direction cosines (\(D.Cos X, D.Cos Y, D.Cos Z\)) and spherical
coordinates (\(\theta, \phi\)) for the
constructs, and one frame with the \(DDISC\) estimates, containing the \(DDISC_{1,2,3}\) parameters for the items
corresponding to the three constructs (1,2,3) (assigned with *c*
in the example above).

`plotD3mirt()`

FunctionThe `plotD3mirt`

function is built on the `rgl`

package (Adler & Murdoch, 2023) for visualization with OpenGL. The
output consists of a three-dimensional interactive RGL device,
displaying vector arrows with the latent dimensions running along the
orthogonal axes centered at zero. If polytomous items are used each item
will have multiple arrows, representing the multiple item step response
functions, running along the exact same direction in the model.

When plotting the `D3mirt`

model with
`plotD3mirt()`

, it is possible to visually observe
statistical violations in the graphical output returned. For instance,
shorter vector arrows indicate weaker discrimination and therefore also
higher amounts of model violations. As another example, if an item
struggles or even fail to describe any of the latent variables in the
model, it can often lead to an extreme stretch of the \(MDIFF\) range. This is comparable to trace
lines turning horizontal in a unidimensional item response theory model.
Some examples of model violations and within-dimensionality will be
given in the illustration below.

To illustrate the utility of `D3mirt`

visualization, a
short example using the `anes0809offwaves`

data set will be
presented. Graphing in default mode by calling `plotd3mirt()`

will return an RGL device that will appear in an external window as a
three-dimensional interactive object that can be rotated manually by the
user. In this illustration, however, all RGL devices are plotted as
inline interactive objects with the help of R markdown. Note, the
`view`

argument is used in all function calls to
`plotD3mirt()`

below. This was only done for optimizing the
size of the RGL device (the default is `view = c(15, 20, 7)`

)
when creating this R Markdown vignette.

```
# Plot RGL device
plotD3mirt(g, view = c(15, 20, 0.6))
```