Table of Contents

# Institutional descrimination: Random-effect of Girl Model

\[ y_{i}^{1} = \Lambda_{i,j}^{1,2} \times \eta_j^{2}+ e_{i}, \]

\[ e_i^{1} \sim N(0, \Theta^{1,1}) \] \[ \eta_{j}^2 \sim N(\alpha^{2}, \Psi^{2,2}) \]

# xxM Model Matrices

## Level-1 (Student): Within matrices

### Residual Covariance Matrix

We have a single dependent variable and hence a single parameter level-1 residual variance \( (\Theta_{1,1}) \) at level-1. Hence, the residual covariance or theta matrix is a (1×1) matrix:

\Theta^{1,1} =

\begin{bmatrix}

\theta^{1,1}_{1,1}

\end{bmatrix}

\]

## Level-2 (School): Within Matrices

At level-2, we have two latent variables: intercept and slope of Girl. Hence, we have two latent means and a covariance matrix.

### Latent Means

With two latent variables, the latent variable mean matrix is a (2×1) matrix:

\alpha^2 =

\begin{bmatrix}

\alpha^2_1 & \\

\alpha^2_2

\end{bmatrix}

\]

\( (\alpha^2_1) \) is the mean of the intercept and \( (\alpha^2_2) \) is the mean of the slope parameter or the average effect of \( (G_{ij}) \) on \( (y_{ij}) \).

### Latent Factor Covariance Matrix

The latent covariance matrix is a (2×2) matrix with two variances and single covariance:

\Psi^{2,2}=

\begin{bmatrix}

\psi_{1,1}^{2,2} & \\

\psi_{2,1}^{2,2} & \psi_{2,2}^{2,2}

\end{bmatrix}

\]

\( (\psi_{1,1}^{2,2}) \) is the variance of the intercept factor representing variability in the intercept of \( (y_{ij}) \) across schools and \( (\psi_{2,2}^{2,2}) \) is the variance of the slope parameter representing between-persons variability in the effect of \( (G_{ij}) \) on \( (y_{ij}) \). Finally, \( (\psi_{2,1}^{2,2}) \) is the covariance between intercept and slope .

## Across level matrices: School to Student

### Across-level factor-loading matrix

We need to capture the effect of level-2 intercept and slope factors on the level-1 dependent variable using a factor-loading matrix with fixed parameters. The factor-loading matrix \( (\Lambda^{1,2}) \) has a single row and two columns (1×2):

\Lambda^{1,2} =

\begin{bmatrix}

1.0_{i,j}^{1,2} & G_{i,j}^{1,2}

\end{bmatrix}

\]

The first column is fixed to 1.0, whereas the second column is fixed to student-specific values of **girl** or \( (G_{ij}) \).

# Path Diagram

(insert diagram)

# Code Listing

## xxM

There is a single **student** level dependent variable **course** and a single independent variable **girl** :

- We do not wish to estimate any factor-loadings. Factor-loadings are fixed to 1.0 and \( (Girl_{ij}) \). Hence, both elements of pattern matrix are zero: \[ \Lambda_{pattern} = \begin{bmatrix} 0 & 0 \end{bmatrix} \]
- We want to fix the first-factor loading to 1.0 (intercept). We use the value matrix to provide the fixed-value of 1.0 for the first factor-loading. The second factor-loading does not have a single fixed value for every observation. Instead each observation \( (i) \) would have its own value for that factor-loading \( (Girl_{ij}) \). Clearly, the value matrix cannot be used for providing individual specific fixed values. Hence, the second element in the value matrix is left as 0.0.
**xxM**ignores it internally. \[ \Lambda_{value} = \begin{bmatrix} 1.0 & 0.0 \end{bmatrix} \] - The job of fixing the factor-loading is left to the label matrix. A label matrix is used to assign labels to each parameter within the matrix. Label matrices can be used to impose equality constraints across matrices. Any two parameters with the same label are constrained to be equal. Label matrices are also used for specifying that a specific parameter is to be fixed to data-values. In this case, the first-label is irrelevant as that parameter has already been fixed to 1.0. We use a descriptive label \( (One) \) as the first label indicating that the loading is fixed to \( 1.0 \). The second factor-loading is the one we are interested in. We want to fix the second factor-loading to the observation specific values of the predictor \( (Girl_{ij}) \). This is accomplished by using a two-part label:
**levelName.predictorName**. In this case, the predictor is a**student**level variable. Hence, the first part of the label is**student**. The second part is the actual predictor name, in this case**girl**. \[ \Lambda_{label} = \begin{bmatrix} One & student.girl \end{bmatrix} \]

### Load xxM

**xxM** library needs to be loaded first.

```
library(xxm)
```

### Prepare R datasets

For this analysis, we use import data from the text file downloaded from XXx. Separate **student** and **school** datasets are created as follows:

```
gcsemv <- read.table(file = "Gcsemv.txt")
# use proper names
names(gcsemv) <- c("school", "student", "girl", "written", "course")
# change type of ID variables to integer
gcsemv$school <- as.integer(gcsemv$school)
gcsemv$student <- 1:nrow(gcsemv)
# Recode (written == -1) and (course == -1 ) as missing
gcsemv$written[gcsemv$written == -1] <- NA
gcsemv$course[gcsemv$course == -1] <- NA
# new dummy variable boy
gcsemv$boy <- 1 - gcsemv$girl
gcsemv$girl <- as.numeric(gcsemv$girl)
gcsemv$boy <- as.numeric(gcsemv$boy)
# school dataset
school <- data.frame(school = unique(gcsemv$school))
# student dataset
student <- gcsemv[c("student", "school", "boy", "girl", "written", "course")]
str(student)
```

```
## 'data.frame': 1905 obs. of 6 variables:
## $ student: int 1 2 3 4 5 6 7 8 9 10 ...
## $ school : int 20920 20920 20920 20920 20920 20920 20920 20920 20920 22520 ...
## $ boy : num 1 0 0 0 1 0 0 1 1 0 ...
## $ girl : num 0 1 1 1 0 1 1 0 0 1 ...
## $ written: num 23 NA 39 36 16 36 49 25 NA 48 ...
## $ course : num NA 71.2 76.8 87.9 44.4 NA 89.8 17.5 32.4 84.2 ...
```

```
str(school)
```

```
## 'data.frame': 73 obs. of 1 variable:
## $ school: int 20920 22520 22710 22738 22908 23208 25241 30474 35270 37224 ...
```

**student** dataset includes the following columns:

**student**Unique student IDs**school**School IDs for**course**dependent varaible**girl**student level predictor

**school** dataset includes the following columns:

**school**Unique school IDs

### Construct R-matrices

For each parameter matrix, construct three related matrices:

**pattern**matrix: A matrix indicating free or fixed parameters.**value**matrix: with start or fixed values for corresponding parameters.**label**matrix: with user friendly label for each parameter.**label**matrix is optional.

```
# student:
theta_11_pat <- matrix(1, 1, 1)
theta_11_val <- diag(2000, 1)
# school:
psi_22_pat <- matrix(c(1, 1, 1, 1), 2, 2)
psi_22_val <- matrix(c(0.1, 0, 0, 0.1), 2, 2)
alpha_2_pat <- matrix(1, 2, 1)
alpha_2_val <- matrix(c(50, -1), 2, 1)
# School-> student lambda_12
lambda_12_pat <- matrix(c(0, 0), 1, 2)
lambda_12_val <- matrix(c(1, 0), 1, 2)
lambda_12_lab <- matrix(c("one", "student.girl"), 1, 2)
```

### Construct model

`xxmModel()`

is used to declare level names. The function returns a model object that is passed as a parameter to subsequent statements. Variable name for the return value can be abything. A

```
xm <- xxmModel(c("student", "school"))
```

```
## A new model was created.
```

### Add submodels

For each declared level `xxmSubmodel()`

is invoked to add corresponding submodel to the model object. The function adds three types of information to the model object:

**parents**declares a list of all parents of the current level.- Level with the independent variable is the
**parent**level. - Level with the dependent variable is the
**child**level.

- Level with the independent variable is the
**variables**declares names of observed dependent (**ys**), observed independent (**xs**) and latent variables (**etas**) for the level.**data**R data object for the current level.

```
### Submodel: Student
xm <- xxmSubmodel(model = xm, level = "student", parents = c("school"), ys = c("course"),
xs = "girl", etas = , data = student)
```

```
## Submodel for level `student` was created.
```

```
### Submodel: School
xm <- xxmSubmodel(model = xm, level = "school", parents = , ys = , xs = , etas = c("inetrcept (boy)",
"girl_minus_boy"), data = school)
```

```
## Submodel for level `school` was created.
```

### Add within-level matrices

For each declared level `xxmWithinMatrix()`

is used to add within-level parameter matrices. For each parameter matrix, the function adds the three matrices constructed earlier:

**pattern****value****label**(optional)

```
### Within-Matrix: Level :: Student, Type :: observed residual-covariance
xm <- xxmWithinMatrix(model = xm, level = "student", type = "theta", pattern = theta_11_pat,
value = theta_11_val)
```

```
##
## 'theta' matrix does not exist and will be added.
## Added `theta` matrix.
```

```
### Within-Matrix: Level :: Student, Type :: latent means
xm <- xxmWithinMatrix(model = xm, level = "school", type = "alpha", pattern = alpha_2_pat,
value = alpha_2_val)
```

```
##
## 'alpha' matrix does not exist and will be added.
## Added `alpha` matrix.
```

```
### Within-Matrix: Level :: School, Type :: latent covariance
xm <- xxmWithinMatrix(model = xm, level = "school", type = "psi", pattern = psi_22_pat,
value = psi_22_val)
```

```
##
## 'psi' matrix does not exist and will be added.
## Added `psi` matrix.
```

### Add across-level matrices

Pairs of levels that share parent-child relationship have regression relationships. `xxmBetweenMatrix()`

is used to add corresponding rergession matrices connecting the two levels.

- Level with the independent variable is the
**parent**level. - Level with the dependent variable is the
**child**level.

For each parameter matrix, the function adds the three matrices constructed earlier:

**pattern****value****label**(optional)

```
### Between-Matrix: Child :: Student, Parent :: School, Type :: factor-loading
xm <- xxmBetweenMatrix(model = xm, parent = "school", child = "student", type = "lambda",
pattern = lambda_12_pat, value = lambda_12_val, label = lambda_12_lab)
```

```
##
## 'lambda' matrix does not exist and will be added.
## Added `lambda` matrix.
```

### Estimate model parameters

Estimation process is initiated by `xxmRun()`

. If all goes well, a q&d summary of the results is printed.

```
xm <- xxmRun(xm)
```

```
## ------------------------------------------------------------------------------
## Estimating model parameters
## ------------------------------------------------------------------------------
## 14719.2840568302
## 14504.3402750169
## 14351.8549969972
## 14195.1600362638
## 14092.7307648035
## 14028.3252514119
## 13997.1833060281
## 13988.3878578752
## 13974.0035334891
## 13972.0476078292
## 13971.6558727026
## 13971.4824057013
## 13971.4456661651
## 13971.4419317272
## 13971.4416183836
## 13971.4416099235
## 13971.4416056594
## 13971.4416054791
## Model converged normally
## nParms: 6
## ------------------------------------------------------------------------------
## *
## 1: student_theta_1_1 :: 169.794 [ 0.000, 0.000]
##
## 2: school_psi_1_1 :: 102.934 [ 0.000, 0.000]
##
## 3: school_psi_1_2 :: -36.527 [ 0.000, 0.000]
##
## 4: school_psi_2_2 :: 47.936 [ 0.000, 0.000]
##
## 5: school_alpha_1_1 :: 69.425 [ 0.000, 0.000]
##
## 6: school_alpha_2_1 :: 7.128 [ 0.000, 0.000]
##
## ------------------------------------------------------------------------------
```

### Estimate profile-likelihood confidence intervals

Once parameters are estimated, confidence intervals are estimated by invoking `xxmCI()`

. Depending on the the number of observations and the complexity of the model, `xxmCI()`

may take a long time to compute. `xxmCI()`

also prints a summary of parameter estimates and CIS.

```
xm <- xxmCI(xm)
```

### View results

A summary of results may be retrived as an R list by a call to `xxmSummary()`

. The returned list has two elements:

**fit**is a list with five elements:**deviance**is \( -2 Log Likelihood \) for the maximum likelihood fit function.**nParameters**is the total number of unique parameters.**nObservations**is the total number of observations across all levels.**aic**is**Akaike’s Information Criterion**or**AIC**computed as \( -2ll + 2*p \).**bic**is**Bayesian Information Criterion**or**BIC**computed as \( -2ll + p*\log(n) \).

**estimates**is a single table of free parameter estimates

All**xxM**parameters have superscripts**{child, parent}**and subscripts**{to, from}**.**xxM**adds a descriptive parameter label if one is not already provided by the user.

```
xxmSummary(xm)
```

```
## $fit
## $fit$deviance
## [1] 13971
##
## $fit$nParameters
## [1] 6
##
## $fit$nObservations
## [1] 1725
##
## $fit$aic
## [1] 13983
##
## $fit$bic
## [1] 14016
##
##
## $estimates
## child parent to from label
## 1 student student course course student_theta_1_1
## 4 school school inetrcept (boy) inetrcept (boy) school_psi_1_1
## 5 school school inetrcept (boy) girl_minus_boy school_psi_1_2
## 6 school school girl_minus_boy girl_minus_boy school_psi_2_2
## 7 school school inetrcept (boy) One school_alpha_1_1
## 8 school school girl_minus_boy One school_alpha_2_1
## estimate low high
## 1 169.794 158.509 182.172
## 4 102.934 68.042 157.854
## 5 -36.527 -73.865 -11.857
## 6 47.936 24.720 86.253
## 7 69.425 66.718 72.100
## 8 7.128 4.888 9.403
```

### Free model object

`xxM`

model object may hog a significant amount of RAM outside of R’s workspace. This memory will automatically be released, when the workspace is cleared by a call to `rm(list=ls())`

or at the end of the R session. Alternatively, it is recommended that `xxmFree()`

may be called to release the memory.

```
xm <- xxmFree(xm)
```

## Proc Mixed

### Import data

```
DATA rdata ;
INFILE "M:\gcseuv.txt"
DSD
LRECL= 27 ;
INPUT
student
school
boy
girl
written
course
;
RUN;
```

For small models such as the current one Proc Mixed is sufficient.

SAS code for a random-slopes model uses a CLASS statement to identify the level-2 units, in this case school. The MODEL statement estimates the fixed-effects \( (\alpha) \). The RANDOM statement specifies that the level-1 intercepts and the effect of level-1 predictor \( (girl_{ij}) \) is allowed to vary across “subject” (i.e., schools). The covariance among the random-effects \( (G) \) is freely estimated (specified by “type = UN”). The \( G \) matrix corresponds to the xxM \( (\Psi) \) matrix. Finally, like all regression models, Proc Mixed estimates the residual variance of the level-1 dependent variable \( (\Theta_{1,1}) \) by default. The important thing to note is that there is one-to-one correspondence between the parameters estimated in Proc Mixed and SEM.

### Proc Mixed

```
PROC MIXED data=rdata method=ML covtest;
CLASS school;
MODEL course = girl/s;
RANDOM int girl/subject=school type=UN G GCORR;
RUN;
```

For large models PROC HPMIXED is required.

### Proc HPMixed

```
PROC HPMIXED data=rdata ;
CLASS school;
MODEL course = girl/s;
RANDOM int girl/subject=school type=UN ;
RUN;
```

```
library(lme4)
```

```
## Loading required package: lattice
## Loading required package: Matrix
```

```
str(student)
```

```
## 'data.frame': 1905 obs. of 6 variables:
## $ student: int 1 2 3 4 5 6 7 8 9 10 ...
## $ school : int 20920 20920 20920 20920 20920 20920 20920 20920 20920 22520 ...
## $ boy : num 1 0 0 0 1 0 0 1 1 0 ...
## $ girl : num 0 1 1 1 0 1 1 0 0 1 ...
## $ written: num 23 NA 39 36 16 36 49 25 NA 48 ...
## $ course : num NA 71.2 76.8 87.9 44.4 NA 89.8 17.5 32.4 84.2 ...
```

```
c.1.lmer.ml <- lmer(course ~ 1 + girl + (1 + girl | school), student, REML = 0)
c.1.lmer.ml
```

```
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: course ~ 1 + girl + (1 + girl | school)
## Data: student
## AIC BIC logLik deviance
## 13983 14016 -6986 13971
## Random effects:
## Groups Name Std.Dev. Corr
## school (Intercept) 10.15
## girl 6.92 -0.52
## Residual 13.03
## Number of obs: 1725, groups: school, 73
## Fixed Effects:
## (Intercept) girl
## 69.43 7.13
```

```
c.1.lmer.reml <- lmer(course ~ 1 + girl + (1 + girl | school), student, REML = 1)
c.1.lmer.reml
```

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: course ~ 1 + girl + (1 + girl | school)
## Data: student
## REML criterion at convergence: 13967
## Random effects:
## Groups Name Std.Dev. Corr
## school (Intercept) 10.24
## girl 7.02 -0.52
## Residual 13.03
## Number of obs: 1725, groups: school, 73
## Fixed Effects:
## (Intercept) girl
## 69.42 7.13
```