Analyzing ILSAs using R

<div class="my-header"><span>Analyzing ILSAs using R</span>
</div>
<div class="my-footer"><span>Francis Huang / huangf@missouri.edu    
</span></div>

---

###Analyzing International Large-scale Assessments using R
#### (w/applied examples)
####Francis L. Huang, PhD

#### 2024-05-06 (updated: 2024-05-13)

---

## NOTE:

There are several excellent R packages that have been developed for analyzing ILSAs:

- https://daniel-caro.com/r-intsvy
- https://ralsa.ineri.org/
- https://www.air.org/project/nces-data-r-project-edsurvey

Focus today is on continuing to expand the concepts explained earlier and applying these concepts using more familiar R workflows (i.e., you are already familiar with R, getting data, running regressions)

---

# Agenda

.bg-washed-yellow.b--gold.ba.bw2.br3.shadow-5.ph4.mt5[
- .green[.f120[Reading in the data]]
- Basic data management
- Getting descriptives
- Fitting regression models
]

---

### Read in the data
- https://oecd.org/pisa/data
- NOTE: I downloaded the **2018** dataset in SPSS format

---

### Load in packages, read in data...

```
library(sjmisc) #for weighted frequencies
library(dplyr) #for data management
library(survey) #for weighted data analysis
dat <- rio::import("c:/data/pisa/CY07_MSU_STU_QQQ.sav")
names(dat) <- tolower(names(dat)) #make lowercase
dat2 <- filter(dat, cnt == "HKG") #respondents from HK only
```

- Once downloaded, load in (install) the necessary packages
- Import the data-- I use the `import` function in the `rio` package (Chan et al., 2021)-- this allows us to retain the SPSS labels in the dataset
- Subset (filter) to Hong Kong data 
- **NOTE** for some reason, the HK school-level data file is completely missing a lot of data. Around 39/40 schools (out of 152) are missing almost *all* data.

NOTE: I the original dataset take a while to download (since it is huge):

I have placed a copy of the HK dataset that can be used in:

```
dat2 <- rio::import("https://github.com/flh3/pubdata/raw/main/mixPV/hkg.sav")
```
.footnote[
Chung-hong Chan, Geoffrey CH Chan, Thomas J. Leeper, and Jason Becker (2021). rio: A Swiss-army knife for data file I/O. R package version 0.5.29.
]

???

sch <- rio::import("C:/Data/pisa/CY07_MSU_SCH_QQQ.sav")
names(sch) <- tolower(names(sch)) #make lowercase
sch2 <- filter(sch, cnt == "HKG") #

sch <- rio::import("C:/Data/pisa/CY07_MSU_SCH_QQQ.sav")
names(sch) <- tolower(names(sch))
sch2 <- filter(sch, cnt == "HKG") #

---

# Agenda

.bg-washed-yellow.b--gold.ba.bw2.br3.shadow-5.ph4.mt5[
- Reading in the data
- .green[.f120[Basic data management]]
- Getting descriptives
- Fitting regression models
]

---

### Find the variables of interest (read the codebook, view the dataset)- select a subset from the entire dataset

```r
sjmisc::find_var(dat2, "Gender")
```

```
  col.nr  var.name                     var.label
1     18 st004d01t Student (Standardized) Gender
```

```r
dat3 <- dplyr::select(dat2, pv1math:pv10math, pv1read:pv10read, 
 gender = st004d01t, escs, cntschid, 
 w_fstuwt, pqschool, st038q03na, immig,
 health = wb150q01ha, stratum)
names(dat3)
```

```
 [1] "pv1math"    "pv2math"    "pv3math"    "pv4math"    "pv5math"    "pv6math"    "pv7math"   
 [8] "pv8math"    "pv9math"    "pv10math"   "pv1read"    "pv2read"    "pv3read"    "pv4read"   
[15] "pv5read"    "pv6read"    "pv7read"    "pv8read"    "pv9read"    "pv10read"   "gender"    
[22] "escs"       "cntschid"   "w_fstuwt"   "pqschool"   "st038q03na" "immig"      "health"    
[29] "stratum"   
```

NOTE: can rename the variable names into something easier to remember! See format above: `newname` = `oldname`.

---

## NOTE on weights:

- The final total student weight is `w_fstuwt`
  - Sum of weights = population size (i.e., 51,101 15 year olds)

```r
sum(dat2$w_fstuwt) #total number of 15-year olds
```

```
[1] 51100.97
```
- Datasets often have what are referred to as replicate weights (in this case, there are 80 replicate weights: `w_fsturwt1` to `w_fsturwt80`)-- this represents another way to analyze the data-- we will not be using those
- There is a lot of documentation that explains this: https://www.oecd.org/pisa/pisaproducts/pisadataanalysismanualspssandsassecondedition.htm 
- Datasets also contain a weight variable referred to as a Senate Weights (`senwt`)-- the sum of the weights in PISA adds up to 5,000 (regardless of educational territory)

```r
sum(dat2$senwt)
```

```
[1] 5000
```

---

#### Need to recode and understand variables more...

```r
sjmisc::frq(dat3, gender)
```

```
Student (Standardized) Gender (gender) <numeric> 
# total N=6037 valid N=6037 mean=1.51 sd=0.50

Value |          Label |    N | Raw % | Valid % | Cum. %
--------------------------------------------------------
    1 |         Female | 2955 | 48.95 |   48.95 |  48.95
    2 |           Male | 3082 | 51.05 |   51.05 | 100.00
    5 |     Valid Skip |    0 |  0.00 |    0.00 | 100.00
    7 | Not Applicable |    0 |  0.00 |    0.00 | 100.00
    8 |        Invalid |    0 |  0.00 |    0.00 | 100.00
    9 |    No Response |    0 |  0.00 |    0.00 | 100.00
 <NA> |           <NA> |    0 |  0.00 |    <NA> |   <NA>
```

```r
dat3$female <- ifelse(dat3$gender == 1, 1, 0)
table(dat3$female)
```

```

0    1 
3082 2955 
```
NOTE: gender is a `numeric` variable. Need to convert. When entered in a regression, text/factors will automatically be dummy coded in R. BUT, can convert to numeric to make things easier (e.g., can run correlations)
---

### ESCS is a continuous measure of socioeconomic status

```r
str(dat3$escs) 
```

```
 num [1:6037] -1.907 -1.977 -2.561 -0.387 -1.231 ...
 - attr(*, "label")= chr "Index of economic, social and cultural status"
 - attr(*, "format.spss")= chr "F8.2"
 - attr(*, "labels")= Named num [1:4] 95 97 98 99
  ..- attr(*, "names")= chr [1:4] "Valid Skip" "Not Applicable" "Invalid" "No Response"
```

```r
hist(dat3$escs)
```

![](ILSA_files/figure-html/unnamed-chunk-7-1.png)

???

#### The `st038` variables are measures of bullying (several are available, showing only one)

```
frq(dat3, st038q03na)
```

---

### A report compares immigrant vs non-immigrant performance in HK

Inspect the variable: NOTE the way we are recoding

```r
frq(dat3, immig)
```

```
Index Immigration status (immig) <numeric> 
# total N=6037 valid N=5814 mean=1.49 sd=0.70

Value |             Label |    N | Raw % | Valid % | Cum. %
-----------------------------------------------------------
    1 |            Native | 3631 | 60.15 |   62.45 |  62.45
    2 | Second-Generation | 1500 | 24.85 |   25.80 |  88.25
    3 |  First-Generation |  683 | 11.31 |   11.75 | 100.00
    5 |        Valid Skip |    0 |  0.00 |    0.00 | 100.00
    7 |    Not Applicable |    0 |  0.00 |    0.00 | 100.00
    8 |           Invalid |    0 |  0.00 |    0.00 | 100.00
    9 |       No Response |    0 |  0.00 |    0.00 | 100.00
 <NA> |              <NA> |  223 |  3.69 |    <NA> |   <NA>
```

```r
#recode non-immigrant vs immigrant / dummy code
dat3$immig2 <- ifelse(dat3$immig > 1, 1, 0)
```

---

### A report compares immigrant vs non-immigrant performance in HK (cont.)

.pull-left[
The following graph shows the difference between immigrant vs nonimmigrant reading outcomes (after controlling for student and school escs)
]

---

### NOTE that even though variables are **labelled**, they are still numeric

```r
frq(dat3, health)
```

```
How is your health? (health) <numeric> 
# total N=6037 valid N=5274 mean=2.01 sd=0.80

Value |          Label |    N | Raw % | Valid % | Cum. %
--------------------------------------------------------
    1 |      Excellent | 1438 | 23.82 |   27.27 |  27.27
    2 |           Good | 2530 | 41.91 |   47.97 |  75.24
    3 |           Fair | 1096 | 18.15 |   20.78 |  96.02
    4 |           Poor |  210 |  3.48 |    3.98 | 100.00
    5 |     Valid Skip |    0 |  0.00 |    0.00 | 100.00
    7 | Not Applicable |    0 |  0.00 |    0.00 | 100.00
    8 |        Invalid |    0 |  0.00 |    0.00 | 100.00
    9 |    No Response |    0 |  0.00 |    0.00 | 100.00
 <NA> |           <NA> |  763 | 12.64 |    <NA> |   <NA>
```

---

### May need to convert certain variables to factors (works with labelled SPSS data)

```r
dat3$health <- rio::factorize(dat3$health)
frq(dat3, health)
```

```
How is your health? (health) <categorical> 
# total N=6037 valid N=5274 mean=2.01 sd=0.80

Value          |    N | Raw % | Valid % | Cum. %
------------------------------------------------
Excellent      | 1438 | 23.82 |   27.27 |  27.27
Good           | 2530 | 41.91 |   47.97 |  75.24
Fair           | 1096 | 18.15 |   20.78 |  96.02
Poor           |  210 |  3.48 |    3.98 | 100.00
Valid Skip     |    0 |  0.00 |    0.00 | 100.00
Not Applicable |    0 |  0.00 |    0.00 | 100.00
Invalid        |    0 |  0.00 |    0.00 | 100.00
No Response    |    0 |  0.00 |    0.00 | 100.00
<NA>           |  763 | 12.64 |    <NA> |   <NA>
```
NOTE: the variable type is now `categorical` (vs `numeric`).

---

### May need to `relevel` factor variables to select a different reference group:

`Poor` is now the reference group.

```r
dat3$health <- relevel(dat3$health, ref = 'Poor') %>% 
  droplevels() #can use droplevels to remove unused levels
table(dat3$health)
```

```

Poor Excellent      Good      Fair 
      210      1438      2530      1096 
```

---

### If you want weighted descriptives, can use the weights in *certain* functions

```r
# unweighted
frq(dat3, health)
```

```
health <categorical> 
# total N=6037 valid N=5274 mean=2.86 sd=0.79

Value     |    N | Raw % | Valid % | Cum. %
-------------------------------------------
Poor      |  210 |  3.48 |    3.98 |   3.98
Excellent | 1438 | 23.82 |   27.27 |  31.25
Good      | 2530 | 41.91 |   47.97 |  79.22
Fair      | 1096 | 18.15 |   20.78 | 100.00
<NA>      |  763 | 12.64 |    <NA> |   <NA>
```
]

```r
# weighted
frq(dat3, health, weights = w_fstuwt)
```

```
health <categorical> 
# total N=44537 valid N=44537 mean=2.85 sd=NA

Value     |     N | Raw % | Valid % | Cum. %
--------------------------------------------
Poor      |  1710 |  3.84 |    3.84 |   3.84
Excellent | 12364 | 27.76 |   27.76 |  31.60
Good      | 21256 | 47.73 |   47.73 |  79.33
Fair      |  9207 | 20.67 |   20.67 | 100.00
<NA>      |     0 |  0.00 |    <NA> |   <NA>
```
]

- At times might not make a big difference

---

### For continuous variables...

```r
mean(dat3$escs, na.rm = TRUE)
```

```
[1] -0.5238875
```

```r
weighted.mean(dat3$escs, dat3$w_fstuwt, na.rm = TRUE)
```

```
[1] -0.5092104
```

---

# Agenda

.bg-washed-yellow.b--gold.ba.bw2.br3.shadow-5.ph4.mt5[
- Reading in the data
- Basic data management
- .green[.f120[Getting descriptives]]
- Fitting regression models
]

---

#### For creating a complete set of weighted descriptives, the `tableone` and `survey` packages are essential

How does `survey` (Lumley, 2010, 2020) work? Need to specify a design first-- this incorporates the design features of the survey.

```r
library(survey)
des <- svydesign(ids = ~cntschid, weights = ~w_fstuwt, data = dat3,
                 stratum = ~strata) 
library(tableone)
t1 <- svyCreateTableOne(vars = c('female', 'st038q03na',
    'health', 'escs', 'immig2'), data = des)
print(t1)
```

```
                        
                         Overall         
  n                      51100.97        
  female (mean (SD))         0.49 (0.50) 
  st038q03na (mean (SD))     1.33 (0.72) 
  health (%)                             
     Poor                  1709.6 ( 3.8) 
     Excellent            12363.6 (27.8) 
     Good                 21256.3 (47.7) 
     Fair                  9207.4 (20.7) 
  escs (mean (SD))          -0.51 (1.04) 
  immig2 (mean (SD))         0.38 (0.49) 
```

### Can also do this as well unweighted (without the survey package)

```r
t0 <- CreateTableOne(vars = c('female', 'st038q03na',
    'health', 'escs', 'immig2'), data = dat3)
print(t0)
```

```
                        
                         Overall      
  n                       6037        
  female (mean (SD))      0.49 (0.50) 
  st038q03na (mean (SD))  1.33 (0.71) 
  health (%)                          
     Poor                  210 ( 4.0) 
     Excellent            1438 (27.3) 
     Good                 2530 (48.0) 
     Fair                 1096 (20.8) 
  escs (mean (SD))       -0.52 (1.02) 
  immig2 (mean (SD))      0.38 (0.48) 
```

---

### If you want to see what is missing...

```r
summary(t0)
```

```

### Summary of continuous variables ###

strata: Overall
              n miss p.miss mean  sd median p25 p75 min max  skew kurt
female     6037    0      0  0.5 0.5    0.0   0 1.0   0   1  0.04 -2.0
st038q03na 6037  443      7  1.3 0.7    1.0   1 1.0   1   4  2.35  4.9
escs       6037  198      3 -0.5 1.0   -0.6  -1 0.2  -7   3 -0.10  0.3
immig2     6037  223      4  0.4 0.5    0.0   0 1.0   0   1  0.51 -1.7

=======================================================================================

### Summary of categorical variables ###

strata: Overall
    var    n miss p.miss     level freq percent cum.percent
 health 6037  763   12.6      Poor  210     4.0         4.0
                         Excellent 1438    27.3        31.2
                              Good 2530    48.0        79.2
                              Fair 1096    20.8       100.0
                                                           
```

---

### Can also do this with correlation tables...

```r
select(dat3, female, escs, 
       st038q03na) %>%
  cor(., use = 'complete.obs') %>%
  round(2)
```

```
           female  escs st038q03na
female       1.00  0.10      -0.15
escs         0.10  1.00      -0.03
st038q03na  -0.15 -0.03       1.00
```
]

```r
library(jtools)
svycor(~female + escs +
    st038q03na, des,
       na.rm = TRUE)
```

```
           female  escs st038q03na
female       1.00  0.10      -0.15
escs         0.10  1.00      -0.03
st038q03na  -0.15 -0.03       1.00
```
NOTE: This requires the `jtools` (Long, 2022) package.
]

In this case, not much of a difference-- but you will not know until you create those tables. That may not always be the case!

.footnote[
Long, J. (2022). jtools: Analysis and Presentation of Social Scientific Data. R package version 2.2.0.
]

---

# Agenda

.bg-washed-yellow.b--gold.ba.bw2.br3.shadow-5.ph4.mt5[
- Reading in the data
- Basic data management
- Getting descriptives
- .green[.f120[Fitting regression models]]
]

---

## Different approaches for running regression models

- Can use single-level models (which can include higher-level predictors) and use cluster robust standard errors (to account for clustering)
- Fit multilevel models
- The regressions must be done for each plausible value and the results combined properly (see figure below; start with imputed data)
- Should use weights
.center[
<img src="https://stefvanbuuren.name/fimd/fig/ch01-miflow-1.png" width="65%" />
]

```r
library(mitools) 
library(mice)
tmp1 <- mitools::withPV(list(maths ~ pv1math + pv2math +
 pv3math + pv4math + pv5math + pv6math + pv7math + 
 pv8math + pv9math + pv10math),
 data = des,
*action = quote(survey::svyglm(maths ~ 1, design = des)),
 rewrite = TRUE)
summary(pool(tmp1)) #from mice package
```

```
         term estimate std.error statistic       df p.value
1 (Intercept) 551.1543  4.759562  115.7994 13708.91       0
```
.center[
<img src="imgs/trend1.jpg" width="70%" />
]

---

#### What about if we do this for reading? (2018)

```r
read1 <- mitools::withPV(list(read ~ pv1read + pv2read +
 pv3read + pv4read + pv5read + pv6read + pv7read + 
 pv8read + pv9read + pv10read),
 data = des,
*action = quote(survey::svyglm(read ~ 1, design = des)),
 rewrite = TRUE)
summary(pool(read1)) #from mice package
```

```
         term estimate std.error statistic     df p.value
1 (Intercept) 524.2831   5.16405  101.5256 298091       0
```
.center[
<img src="imgs/trend1.jpg" width="70%" />
]

---

```r
tmp2 <- mitools::withPV(list(maths ~ pv1math + pv2math +
 pv3math + pv4math + pv5math + pv6math + pv7math + 
 pv8math + pv9math + pv10math), data = des,
*action = quote(survey::svyglm(maths ~ female, design = des)),
 rewrite = TRUE)
summary(pool(tmp2))
```

```
         term   estimate std.error  statistic       df   p.value
1 (Intercept) 548.459882  5.388303 101.787131 140.2126 0.0000000
2      female   5.537028  5.009030   1.105409 122.7376 0.2711444
```

---

#### A few slides ago, showed the reading performance gap of immigrant vs nonimmigrant students after controlling to student and school SES (9 point difference)

- Need to first create a group-level aggregate of SES at the school level:

```r
library(MLMusingR)
dat3$escs_sch <- group_mean(dat3$escs, dat3$cntschid)
```

- Need to update the `des` object since have added a new variable to the dataset
- `escs_sch` is referred to as a school-level variable-- this is different from the student-level variable
  - The same value repeats per school

---

### Regression results, including immigration status and both student- and school-level escs as predictors

```r
des <- svydesign(ids = ~cntschid, weights = ~w_fstuwt, data = dat3) 
tmp3 <- mitools::withPV(list(read ~ pv1read + pv2read +
 pv3read + pv4read + pv5read + pv6read + pv7read + 
 pv8read + pv9read + pv10read), data = des,
*action = quote(survey::svyglm(read ~ immig2 + escs + escs_sch, design = des)),
 rewrite = TRUE)
summary(pool(tmp3))
```

```
         term   estimate std.error  statistic       df      p.value
1 (Intercept) 553.373765  5.408380 102.317840 144.3773 0.000000e+00
2      immig2   8.823202  4.483261   1.968032 136.3429 5.109335e-02
3        escs   3.525814  1.560557   2.259330 108.2792 2.586512e-02
4    escs_sch  55.813817  7.019104   7.951701 142.8540 5.082601e-13
```
- From the report:
> "There is no statistically significant difference in reading performance between immigrant and nonimmigrant students in Hong Kong (China). After accounting for students' and schools' socio-economic profile the difference shrank to nine score points." (p. 6)
---

### Another way: 
1. Create a tall dataset from the wide dataset: one outcome per row (i.e., read)
2. Split the data per pv (into a list)
3. Run a regression per item in a list
4. Pool the results

```r
library(tidyr) #for reshaping
tall <- pivot_longer(dat3, pv1read:pv10read, names_to = 'pv', values_to = 'read')
listofdata <- split(tall, tall$pv)
library(estimatr) #for lm_robust
listofres <- lapply(listofdata, function(x)
  lm_robust(read ~ immig2 + escs + escs_sch, data = x,
  weights = w_fstuwt, cluster = cntschid)) #can use se_type = 'stata' too
summary(pool(listofres))
```

```
         term   estimate std.error  statistic       df      p.value
1 (Intercept) 553.373765  5.491749 100.764582 35.19116 0.000000e+00
2      immig2   8.823202  4.571647   1.929983 33.92185 6.200541e-02
3        escs   3.525814  1.564604   2.253487 30.05278 3.168059e-02
4    escs_sch  55.813817  7.187124   7.765807 34.92349 4.097163e-09
```
---

#### Fitting the models per PV shows differences... (showing pv 10, 1 to 5)

<table style="NAborder-bottom: 0; width: auto !important; margin-left: auto; margin-right: auto;" class="table">
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:center;"> pv10read </th>
   <th style="text-align:center;"> pv1read </th>
   <th style="text-align:center;"> pv2read </th>
   <th style="text-align:center;"> pv3read </th>
   <th style="text-align:center;"> pv4read </th>
   <th style="text-align:center;"> pv5read </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:center;"> 553.208*** </td>
   <td style="text-align:center;"> 553.113*** </td>
   <td style="text-align:center;"> 553.639*** </td>
   <td style="text-align:center;"> 553.797*** </td>
   <td style="text-align:center;"> 552.998*** </td>
   <td style="text-align:center;"> 553.573*** </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (5.323) </td>
   <td style="text-align:center;"> (5.585) </td>
   <td style="text-align:center;"> (5.674) </td>
   <td style="text-align:center;"> (5.424) </td>
   <td style="text-align:center;"> (5.467) </td>
   <td style="text-align:center;"> (5.439) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> immig2 </td>
   <td style="text-align:center;"> 9.564* </td>
   <td style="text-align:center;"> 9.591* </td>
   <td style="text-align:center;"> 9.209* </td>
   <td style="text-align:center;"> 7.674+ </td>
   <td style="text-align:center;"> 9.635* </td>
   <td style="text-align:center;"> 7.760+ </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (4.493) </td>
   <td style="text-align:center;"> (4.536) </td>
   <td style="text-align:center;"> (4.408) </td>
   <td style="text-align:center;"> (4.589) </td>
   <td style="text-align:center;"> (4.510) </td>
   <td style="text-align:center;"> (4.467) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> escs </td>
   <td style="text-align:center;"> 3.241* </td>
   <td style="text-align:center;"> 3.096* </td>
   <td style="text-align:center;"> 3.793* </td>
   <td style="text-align:center;"> 3.260* </td>
   <td style="text-align:center;"> 3.544* </td>
   <td style="text-align:center;"> 4.531** </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (1.493) </td>
   <td style="text-align:center;"> (1.453) </td>
   <td style="text-align:center;"> (1.487) </td>
   <td style="text-align:center;"> (1.575) </td>
   <td style="text-align:center;"> (1.424) </td>
   <td style="text-align:center;"> (1.476) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> escs_sch </td>
   <td style="text-align:center;"> 57.171*** </td>
   <td style="text-align:center;"> 56.167*** </td>
   <td style="text-align:center;"> 55.768*** </td>
   <td style="text-align:center;"> 55.505*** </td>
   <td style="text-align:center;"> 55.776*** </td>
   <td style="text-align:center;"> 54.072*** </td>
  </tr>
  <tr>
   <td style="text-align:left;box-shadow: 0px 1px">  </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (6.862) </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (7.345) </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (7.337) </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (7.099) </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (7.095) </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (7.141) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> R2 </td>
   <td style="text-align:center;"> 0.135 </td>
   <td style="text-align:center;"> 0.128 </td>
   <td style="text-align:center;"> 0.130 </td>
   <td style="text-align:center;"> 0.128 </td>
   <td style="text-align:center;"> 0.131 </td>
   <td style="text-align:center;"> 0.128 </td>
  </tr>
</tbody>
<tfoot><tr><td style="padding: 0; " colspan="100%">
<sup></sup> + p &lt; 0.1, * p &lt; 0.05, ** p &lt; 0.01, *** p &lt; 0.001</td></tr></tfoot>
</table>

---

### Showing results by PV shows...

- Conclusions can differ by PV
  - In the analysis per PV, some were statistically significant-- depending on the PV chosen!

- p-values in the pooled results > .05 (i.e., .051, .062)- careful about interpreting p-values too much

---

### Note, in all the models we have fit, the student-level coefficients are a mix of effects

- There is a portion due to the student
- There is a portion due to the school attended
- We can disentangle that but I did not do that to keep results consistent with the reports that have been produced
  - One way of doing that is to group-center level-1 predictors
  - Another way is to include the group mean as a predictor at level 2
  - This applies to both single- and multilevel models
- Group centering is merely subtracting the group mean from the original variable- can be done with continuous and binary variables
- Sometimes referred to as centering within context (CWC)

---

### Model math achievement both ways using GLM and MLM using group-centered variables

```r
library(MLMusingR) #has the group_center function
math.glm <- mitools::withPV(list(maths ~ pv1math + pv2math +
 pv3math + pv4math + pv5math + pv6math + pv7math + 
 pv8math + pv9math + pv10math), data = des,
 action = quote(survey::svyglm(maths ~ group_center(female, cntschid), 
 design = des)), 
 rewrite = TRUE)

math.mlm <- mixPV(pv1math + pv2math +
 pv3math + pv4math + pv5math + pv6math + pv7math + 
 pv8math + pv9math + pv10math ~ group_center(female, cntschid) +
 (1|cntschid), silent = TRUE,
 data = dat3, weights = c('w_fstuwt', 'one')) 
```

---

### Comparison of results (uncentered, centered GLM, centered MLM)

---

### Are weights needed?

- If the use of weights were 'cost-free' (i.e., you do not give up anything using weights), then the use of weights would be routine
- Some disciplines do not generally tend to use weights (e.g., econometrics)
- However, the use of weights does make sense as it 'respects' the study design
- Sometimes, weights may not be *informative* in the sense that results do not change
- However, in those cases, the variance of the coefficients become larger (inflated) unnecessarily (i.e., results in less power)
- Some tests have been formulated to test if weights are informative (see Bollen et al., 2016 for a very good and readable overview)

.footnote[
Bollen, K. A., Biemer, P. P., Karr, A. F., Tueller, S., & Berzofsky, M. E. (2016). Are survey weights needed? A review of diagnostic tests in regression analysis. Annual Review of Statistics and Its Application, 3(1), 375–392. https://doi.org/10.1146/annurev-statistics-011516-012958
]

---

#### There are some basic tests and are simple to understand (in this case, the weights are informative)

- There are other functions as well that are helpful (using one pv as an example)

```r
library(jtools) #fit model, use weights_test
tst <- lm(pv1math ~ female + escs, data = dat3) #model to test
weights_tests(model = tst, weights = w_fstuwt, data = dat3) #run the test
```

```
DuMouchel-Duncan test of model change with weights

F(3,5833) = 47.123
p = 0

Lower p values indicate greater influence of the weights.

Standard errors: OLS
-----------------------------------------------------
                          Est.   S.E.   t val.      p
--------------------- -------- ------ -------- ------
(Intercept)             616.63   6.56    94.04   0.00
female                  -31.91   9.18    -3.48   0.00
w_fstuwt                 -5.77   0.74    -7.84   0.00
escs                     -2.14   4.13    -0.52   0.61
female:w_fstuwt           3.20   1.05     3.05   0.00
w_fstuwt:escs             2.55   0.46     5.61   0.00
-----------------------------------------------------

---
Pfeffermann-Sverchkov test of sample weight ignorability

Residual correlation = -0.13, p = 0.00
Squared residual correlation = 0.08, p = 0.00
Cubed residual correlation = -0.11, p = 0.00

A significant correlation may indicate biased estimates
in the unweighted model.
```

---

### One test (DuMouchel-Duncan, 1983) fits the basic model and then compares this to a model with weights interacting with the predictors...

- If the model fits better with the interactions, the weights are informative and should be used
- See prior slide-- results are the same F(3, 5833) = 47.12, p < .001

```r
tst2 <- update(tst, . ~ . * w_fstuwt)
anova(tst, tst2)
```

```
Analysis of Variance Table

Model 1: pv1math ~ female + escs
Model 2: pv1math ~ female + escs + w_fstuwt + female:w_fstuwt + escs:w_fstuwt
  Res.Df      RSS Df Sum of Sq      F    Pr(>F)    
1   5836 47526832                                  
2   5833 46402224  3   1124608 47.123 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

- The Pfeffermann-Sverchkov test uses the residuals (at the bottom of the other page, but cannot be seen on the slide)

---

## Summary

- There are different approaches for using weighted data (GLM, MLM)
- Different packages can handle the use the plausible values
- Be mindful of the types of weights available
- Tests can be performed in order to test the need for the use of weights