Model with Binomially Distributed Responses
We use a simulated dataset from the PLmixed package for this example. The observations are binomial responses representing an ability measurement. The first few lines are shown below.
library(PLmixed)
head(IRTsim)
#> sid school item y
#> 1.1 1 1 1 1
#> 1.2 1 1 2 1
#> 1.3 1 1 3 1
#> 1.4 1 1 4 0
#> 1.5 1 1 5 1
#> 2.1 2 1 1 1Each student is identified by a student id sid, and each
school with a school id given by the school variable. For
each student, five item measurements have been made. We assume that the
student’s performance depends both on the students ability as well as on
the school the student attends. Having the outline of GALAMMs from the
introductory
vignette in mind, we assume a binomial response model with a logit
link, yielding for the \(i\)th
measurement of the \(j\)th student in
the \(k\)th school,
\[ \text{P}(y_{ijk} = 1 | \mathbf{x}_{ijk}, \boldsymbol{\eta}_{jk}) = \frac{\exp(\nu_{ijk})}{1 + \exp(\nu_{ijk})} \]
The latent variable vector is given by \(\boldsymbol{\eta}_{jk} = (\eta_{j}, \eta_{jk})^{T}\), where the last element is the effect of the school which the \(j\)th student attends. The nonlinear predictor is given by
\[ \nu_{ijk} = \mathbf{x}_{ijk}^{T} \boldsymbol{\beta} + \mathbf{x}_{ijk}^{T}\boldsymbol{\lambda} (\eta_{j} + \eta_{jk}) \]
where \(\mathbf{x}_{ij}\) is a vector containing dummy variables for items, \(\eta_{j}\) is a latent variable describing the underlying ability of student \(j\), and \(\eta_{k}\) is a latent variable describing the effect of school \(k\). The factor loading \(\lambda\) describes how the sum of the two latent variables impacts the nonlinear predictor, and hence the probability of correct response. We note that using a common \(\lambda\) for two latent variables like this is not necessarily a good model, but it makes this introductory example easier to follow. We will extend the model later.
The structural model is simply
\[ \boldsymbol{\eta}_{jk} = \boldsymbol{\zeta}_{jk} \sim N_{2}(\mathbf{0}, \boldsymbol{\Psi}), \]
where \(N(\mathbf{0}, \boldsymbol{\Psi})\) denotes a bivariate normal distribution with mean zero covariance matrix \(\boldsymbol{\Psi}\). The covariance matrix is assumed to be diagonal.
We confirm that item has five levels. This means that
\(\boldsymbol{\lambda}\) is a vector
with five elements.
For identifiability, we fix the first element of \(\boldsymbol{\lambda}\) to one. The rest
will be freely estimated. We do this by defining the following matrix.
Any numeric value implies that the element is fixed, whereas
NA implies that the element is unknown, and will be
estimated.
(loading_matrix <- matrix(c(1, NA, NA, NA, NA), ncol = 1))
#> [,1]
#> [1,] 1
#> [2,] NA
#> [3,] NA
#> [4,] NA
#> [5,] NAWe fit the model as follows:
mod <- galamm(
formula = y ~ item + (0 + ability | school / sid),
data = IRTsim,
family = binomial,
load.var = "item",
factor = "ability",
lambda = loading_matrix
)A couple of things in the model formula are worth pointing out. First
the part (0 + ability | school / sid) in the model formula
specifies that student ability varies between students within schools.
It corresponds to the term \(\mathbf{x}_{ijk}^{T}\boldsymbol{\lambda} (\eta_{j}
+ \eta_{jk})\) in the mathematical specification of the model.
The variable ability is not part of the IRTsim
dataframe, but is instead specified in the argument
factor = "ability". The argument
load.var = "item" specifies that all rows in the dataframe
with the same value of “item” should get the same element of \(\boldsymbol{\lambda}\), and hence it
defines the dummy variable \(\mathbf{x}_{ijk}\). Finally,
lambda = loading_matrix provides the matrix of factor
loadings. Note that we must explicitly add a zero in
(0 + ability | school / sid) to avoid having a random
intercept estimated in addition to the effect for each value of “item”;
such a model would not be identified. The fixed effect part of the model
formula, which is simply item, specifies the term \(\mathbf{x}_{ijk}^{T}
\boldsymbol{\beta}\).
We can start by inspecting the fitted model:
summary(mod)
#> GALAMM fit by maximum marginal likelihood.
#> Formula: y ~ item + (0 + ability | school/sid)
#> Data: IRTsim
#>
#> AIC BIC logLik deviance df.resid
#> 2966.4 3030.5 -1472.2 2372.3 2489
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.8884 -0.7335 0.4266 0.6235 3.0410
#>
#> Lambda:
#> ability SE
#> lambda1 1.0000 .
#> lambda2 0.7370 0.1456
#> lambda3 0.9351 0.1872
#> lambda4 0.6069 0.1261
#> lambda5 0.5860 0.1163
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> sid:school ability 1.466 1.211
#> school ability 1.300 1.140
#> Number of obs: 2500, groups: sid:school, 500; school, 26
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.5112 0.2617 1.953 0.050777
#> item2 0.3256 0.1791 1.818 0.069068
#> item3 -0.4491 0.1623 -2.768 0.005644
#> item4 0.4930 0.1924 2.562 0.010402
#> item5 0.4585 0.1922 2.385 0.017067The fixef method lets us consider the fixed effects:
fixef(mod)
#> (Intercept) item2 item3 item4 item5
#> 0.5112006 0.3255941 -0.4491070 0.4929827 0.4584897We can also get Wald type confidence intervals for the fixed effects.
confint(mod, parm = "beta")
#> 2.5 % 97.5 %
#> (Intercept) -0.001727334 1.0241286
#> item2 -0.025430062 0.6766183
#> item3 -0.767138526 -0.1310755
#> item4 0.115869898 0.8700955
#> item5 0.081749228 0.8352303We can similarly extract the factor loadings.
factor_loadings(mod)
#> ability SE
#> lambda1 1.0000000 NA
#> lambda2 0.7370254 0.1455782
#> lambda3 0.9351105 0.1872049
#> lambda4 0.6069065 0.1260943
#> lambda5 0.5859914 0.1162977And we can find confidence intervals for them. Currently, only Wald type confidence intervals are available. Be aware that such intervals may have poor coverage properties.
confint(mod, parm = "lambda")
#> 2.5 % 97.5 %
#> lambda1 0.4516974 1.0223534
#> lambda2 0.5681956 1.3020254
#> lambda3 0.3597662 0.8540468
#> lambda4 0.3580520 0.8139307We can also show a diagnostic plot, although for a binomial model like this it is not very informative.
