Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions vignettes-raw/glmm_factor.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -41,36 +41,36 @@ head(IRTsim)
Each 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](https://docs.ropensci.org/galamm/articles/galamm.html) 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})}
\text{P}(y_{ijk} = 1 | \mathbf{x}_{ijk}, \pmb{\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
The latent variable vector is given by $\pmb{\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})
\nu_{ijk} = \mathbf{x}_{ijk}^{T} \pmb{\beta} + \mathbf{x}_{ijk}^{T}\pmb{\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}),
\pmb{\eta}_{jk} = \pmb{\zeta}_{jk} \sim N_{2}(\mathbf{0}, \pmb{\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.
where $N(\mathbf{0}, \pmb{\Psi})$ denotes a bivariate normal distribution with mean zero covariance matrix $\pmb{\Psi}$. The covariance matrix is assumed to be diagonal.

```{r}
IRTsim$item <- factor(IRTsim$item)
```

We confirm that `item` has five levels. This means that $\boldsymbol{\lambda}$ is a vector with five elements.
We confirm that `item` has five levels. This means that $\pmb{\lambda}$ is a vector with five elements.

```{r}
table(IRTsim$item)
```

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.
For identifiability, we fix the first element of $\pmb{\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.

```{r}
(loading_matrix <- matrix(c(1, NA, NA, NA, NA), ncol = 1))
Expand All @@ -89,7 +89,7 @@ mod <- galamm(
)
```

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}$.
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}\pmb{\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 $\pmb{\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} \pmb{\beta}$.

We can start by inspecting the fitted model:

Expand Down
10 changes: 5 additions & 5 deletions vignettes-raw/latent_observed_interaction.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ y_{2} \\
y_{3}
\end{pmatrix}
=
\boldsymbol{\beta}_{0} +
\pmb{\beta}_{0} +
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
Expand All @@ -65,10 +65,10 @@ y_{3}
0 \\
x \beta
\end{pmatrix}
+ \boldsymbol{\epsilon}.
+ \pmb{\epsilon}.
$$

In this equation $\boldsymbol{\beta}_{0} \in \mathbb{R}^{3}$ is a vector of intercepts, $\eta$ is a latent variable, the loading of the latent variable onto the first measurement $y_{1}$ is fixed to 1 for identifiability, $\lambda_{2}$ is the loading of the latent variable onto the second measurement $y_{2}$, $\lambda_{3}$ is the main effect of the latent variable on the response of interest $y_{3}$, $\beta$ is the effect of the observed covariate $x$ on $y_{3}$, and $\lambda_{4}$ is the interaction effect of $x$ and $\eta$ on $y_{3}$. We assume that the residuals $\boldsymbol{\epsilon}$ are independently and identically normally distributed; this assumption is valid in this simulated case, but note that since the response $y_{3}$ is qualitatively different from the measurements $y_{1}$ and $y_{2}$, this assumption will in general not hold, and a [heteroscedastic measurement model](https://docs.ropensci.org/galamm/articles/lmm_heteroscedastic.html) should be used, or a [model with mixed response types](https://docs.ropensci.org/galamm/articles/mixed_response.html). For a more detailed explanation of this way of formulating latent variable models in matrix form we refer to the first four pages of @rockwoodEstimatingComplexMeasurement2019.
In this equation $\pmb{\beta}_{0} \in \mathbb{R}^{3}$ is a vector of intercepts, $\eta$ is a latent variable, the loading of the latent variable onto the first measurement $y_{1}$ is fixed to 1 for identifiability, $\lambda_{2}$ is the loading of the latent variable onto the second measurement $y_{2}$, $\lambda_{3}$ is the main effect of the latent variable on the response of interest $y_{3}$, $\beta$ is the effect of the observed covariate $x$ on $y_{3}$, and $\lambda_{4}$ is the interaction effect of $x$ and $\eta$ on $y_{3}$. We assume that the residuals $\pmb{\epsilon}$ are independently and identically normally distributed; this assumption is valid in this simulated case, but note that since the response $y_{3}$ is qualitatively different from the measurements $y_{1}$ and $y_{2}$, this assumption will in general not hold, and a [heteroscedastic measurement model](https://docs.ropensci.org/galamm/articles/lmm_heteroscedastic.html) should be used, or a [model with mixed response types](https://docs.ropensci.org/galamm/articles/mixed_response.html). For a more detailed explanation of this way of formulating latent variable models in matrix form we refer to the first four pages of @rockwoodEstimatingComplexMeasurement2019.

The structural model is simply $\eta = \zeta \sim N(0, \psi)$, where $\psi$ is its variance.

Expand Down Expand Up @@ -106,7 +106,7 @@ y_{2} \\
y_{3}
\end{pmatrix}
=
\boldsymbol{\beta}_{0} +
\pmb{\beta}_{0} +
\begin{pmatrix}
1 \\
\lambda_{2} \\
Expand All @@ -119,7 +119,7 @@ y_{3}
0 \\
x \beta
\end{pmatrix}
+ \boldsymbol{\epsilon}.
+ \pmb{\epsilon}.
$$

This way of writing shows more explicitly which factor loadings are connected with which observation. In order to fit this model with `galamm`, we must provide formulas for the terms in the loading matrix
Expand Down
14 changes: 7 additions & 7 deletions vignettes-raw/lmm_factor.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -51,10 +51,10 @@ levels(KYPSsim$time)
Since students attending a given middle school not necessarily attend the same high school, we have a model with crossed random effects. We use the model of @jeonProfileLikelihoodApproachEstimating2012, whose measurement part can be formulated as

$$
y_{tsmh} = \beta_{0} + \sum_{t'=2}^{4} d_{tt'}\beta_{t'} + \mathbf{d}_{t}^{T} \left(\boldsymbol{\lambda}_{m} \eta_{m} + \boldsymbol{\lambda}_{h} \eta_{h}\right) + \eta_{s} + \epsilon_{tsmh},
y_{tsmh} = \beta_{0} + \sum_{t'=2}^{4} d_{tt'}\beta_{t'} + \mathbf{d}_{t}^{T} \left(\pmb{\lambda}_{m} \eta_{m} + \pmb{\lambda}_{h} \eta_{h}\right) + \eta_{s} + \epsilon_{tsmh},
$$

where $\mathbf{d}_{t} = (d_{t1},d_{t2},d_{t3},d_{t4})^{T}$ is a vector whose $t$th element equals one and all other elements equal zero. $\beta_{0}$ is an intercept, and $\beta_{2}$, $\beta_{3}$, and $\beta_{4}$ are the effects of timepoints 2, 3, and 4. $\eta_{m}$ and $\eta_{h}$ are the "teacher effects" of middle school $m$ and high school $s$, respectively, $\eta_{s}$ is the latent level for student $s$, and $\epsilon_{tsmh}$ is a residual term. $\boldsymbol{\lambda}_{m}$ and $\boldsymbol{\lambda}_{h}$ are factor loadings (called "persistence parameters" by @mccaffreyModelsValueAddedModeling2004 and @jeonProfileLikelihoodApproachEstimating2012) specifying how the teacher effects for middle school and high school impact the self esteem measurement. Since students attend high school after middle school, the measurements of self esteem in middle school are assumed not to be affected by high school, and hence the first two elements of $\boldsymbol{\lambda}_{h}$ are set to zero. The first nonzero element is set to zero for identifiability, so we have $\boldsymbol{\lambda}_{h} = (0, 0, 1, \lambda_{h4})^{T}$. Conversely, we allow for middle school to have an effect of measurements in high school, so $\boldsymbol{\lambda}_{m} = (1, \lambda_{m2}, \lambda_{m3}, \lambda_{m4})^{T}$, with the first element set to zero for identifiability. The residuals are assumed normally distributed, $\epsilon_{tsmh} \sim N(0, \phi)$.
where $\mathbf{d}_{t} = (d_{t1},d_{t2},d_{t3},d_{t4})^{T}$ is a vector whose $t$th element equals one and all other elements equal zero. $\beta_{0}$ is an intercept, and $\beta_{2}$, $\beta_{3}$, and $\beta_{4}$ are the effects of timepoints 2, 3, and 4. $\eta_{m}$ and $\eta_{h}$ are the "teacher effects" of middle school $m$ and high school $s$, respectively, $\eta_{s}$ is the latent level for student $s$, and $\epsilon_{tsmh}$ is a residual term. $\pmb{\lambda}_{m}$ and $\pmb{\lambda}_{h}$ are factor loadings (called "persistence parameters" by @mccaffreyModelsValueAddedModeling2004 and @jeonProfileLikelihoodApproachEstimating2012) specifying how the teacher effects for middle school and high school impact the self esteem measurement. Since students attend high school after middle school, the measurements of self esteem in middle school are assumed not to be affected by high school, and hence the first two elements of $\pmb{\lambda}_{h}$ are set to zero. The first nonzero element is set to zero for identifiability, so we have $\pmb{\lambda}_{h} = (0, 0, 1, \lambda_{h4})^{T}$. Conversely, we allow for middle school to have an effect of measurements in high school, so $\pmb{\lambda}_{m} = (1, \lambda_{m2}, \lambda_{m3}, \lambda_{m4})^{T}$, with the first element set to zero for identifiability. The residuals are assumed normally distributed, $\epsilon_{tsmh} \sim N(0, \phi)$.

Written out for each of the four timepoints, the model becomes

Expand Down Expand Up @@ -93,7 +93,7 @@ $$

where $N_{3}(a, b)$ denotes a trivariate normal distribution with mean $a$ and covariance $b$.

In order to fit the model with galamm, we use syntax similar to PLmixed, and start by defining the loading matrix. The first column contains $\boldsymbol{\lambda}_{m}$ and the second column contains $\boldsymbol{\lambda}_{h}$. Numerical values in this matrix means that the entry is fixed to the given value, whereas `NA` means that the value is unknown, and should be estimated.
In order to fit the model with galamm, we use syntax similar to PLmixed, and start by defining the loading matrix. The first column contains $\pmb{\lambda}_{m}$ and the second column contains $\pmb{\lambda}_{h}$. Numerical values in this matrix means that the entry is fixed to the given value, whereas `NA` means that the value is unknown, and should be estimated.


```{r}
Expand Down Expand Up @@ -170,7 +170,7 @@ A diagnostic plot of residuals versus predicted values also looks acceptable, al
plot(mod, abline = c(0, 0))
```

We can also compare the estimated model to a model with constrained factor loadings. In particular, we could assume that the teacher effect in middle school has no effect on self esteem measured during high school, by setting the last two elements of $\boldsymbol{\lambda}_{m}$ to zero. We would then have the following loading matrix.
We can also compare the estimated model to a model with constrained factor loadings. In particular, we could assume that the teacher effect in middle school has no effect on self esteem measured during high school, by setting the last two elements of $\pmb{\lambda}_{m}$ to zero. We would then have the following loading matrix.

```{r}
(loading_matrix_constr1 <- rbind(
Expand Down Expand Up @@ -303,7 +303,7 @@ y_{12tsc} \\
\eta_{6s}^{(s)} \\
\eta_{7c}^{(c)} \\
\end{pmatrix} +
\boldsymbol{\epsilon}_{tsc}
\pmb{\epsilon}_{tsc}
$$

In brief, $\eta_{1t}^{(t)}$ and $\eta_{2t}^{(t)}$ are the teacher effects, $\eta_{3s}^{(s)}$ and $\eta_{4s}^{(s)}$ are teacher's perception of the students on the trait, $\eta_{5s}^{(s)}$ and $\eta_{6s}^{(s)}$ are the students' perception of themselves on the trait, and $\eta_{7c}^{(c)}$ is the classroom effect. The factor loadings are the "regression coefficients" for regressing the observed items onto these latent traits. The subscripts $t$, $s$, and $c$ indicate teacher, student, and class, respectively.
Expand Down Expand Up @@ -339,7 +339,7 @@ $$
\zeta_{1t}^{(t)} \\
\zeta_{2t}^{(t)} \\
\end{pmatrix}
\sim N_{2}(\mathbf{0}, \boldsymbol{\Psi}^{(t)}),
\sim N_{2}(\mathbf{0}, \pmb{\Psi}^{(t)}),
$$

$$
Expand All @@ -349,7 +349,7 @@ $$
\zeta_{5s}^{(s)} \\
\zeta_{6s}^{(s)} \\
\end{pmatrix}
\sim N_{4}(\mathbf{0}, \boldsymbol{\Psi}^{(s)}),
\sim N_{4}(\mathbf{0}, \pmb{\Psi}^{(s)}),
$$

$$
Expand Down
2 changes: 1 addition & 1 deletion vignettes-raw/lmm_heteroscedastic.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ library(galamm)
At the moment, galamm supports group-wise heteroscedasticity in Gaussian response models. Referring to the model formulation outlined in the [introductory vignette](https://docs.ropensci.org/galamm/articles/galamm.html), the response model and nonlinear predictor can be easily combined in this particular case, to give

$$
y_{i} = \sum_{s=1}^{S} f_{s}\left(\mathbf{x}_{i}\right) + \sum_{l=2}^{L}\sum_{m=1}^{M_{l}} \eta_{m}^{(l)} \mathbf{z}^{(l)}_{im}{}^{'}\boldsymbol{\lambda}_{m}^{(l)} + \epsilon_{g(i)},
y_{i} = \sum_{s=1}^{S} f_{s}\left(\mathbf{x}_{i}\right) + \sum_{l=2}^{L}\sum_{m=1}^{M_{l}} \eta_{m}^{(l)} \mathbf{z}^{(l)}_{im}{}^{'}\pmb{\lambda}_{m}^{(l)} + \epsilon_{g(i)},
$$

where subscript $i$ refers to the $i$th elementary unit of observation, i.e., the $i$th row in the dataframe. $g(i)$ refers to the group to which the $i$th observation belongs, with each grouping having a separately estimated residual variance, $\epsilon_{g} \sim N(0, \sigma_{g}^{2})$.
Expand Down
Loading
Loading