---
title: "Edps 587: Mixed Effects Location-Scale Models for Nels Data"
author: "Carolyn J. Anderson"
date: '2023-03-17'
output: 
   html_document:
     toc: yes
---
Here fit location-scale model to 23 schools from nels.  The R belows will run the program mixreglsb.exe from R.  The code makes use of the function R_mixregls.txt which is a function that I wrote to help you set up the input to mixreglsb.exe and than run the program.

The mixreglsb.txt use software described in the following paper.  The R described in the paper no longer works; however, you can download programs and examples from the paper's web-site.  There are quite a few number of files, but you only need mixreglsb.exe.

Hedeker, D, & Nordgren, R. (2013), MIXREGLS: A Program for Mixed-Effects Location Scale Analysis. Journal of Statistical Software, 52. http://www.jstatsoft.org/

I suggest that you read the above paper to better understand what is going on.

The mixreglsb.exe only fits random intercept models, but an extended verision is described in the paper Nordrend, Hedeker, Donton & Yang (2019) in Statistics in Medicine.  There is an appendix that contians information for fitting the extended version in SAS/NLMIXED and merlin in STATA.  A R version is in development at the time the paper was written.

A few things if you run the R_mixregls function:

1.    Data are numeric.  If you have any varibles with character values these should be dummy, effect or otherwise coded.
2.    The program mixreglsb.exe runs on a PC.
3.    I have not written in what to do with missing data; however, if you
      have missing data you could create the .def file and then edit to 
      include missing values codes.
4.    You should include the formula.tools and stringr libaries


```{r libraries, message=FALSE, warnings=FALSE}
library(formula.tools)
library(stringr)
```

This is where my data and the mixreglsb.exe files live:
```{r setworking_directory}
setwd("D:/Dropbox/edps587/lectures/8 modelbuilding/MIXREGLS/school23")
```

I will use as source the function
```{r source_of_R_function}
source("D:/Dropbox/edps587/lectures/8 modelbuilding/MIXREGLS/school23/R_mixregls.R")
```

# Syntax

Most of the options are set to defaults but can easily be changed as will be illustrated below.  Also all of these are described in Hekeke & Nordgren (2013)

```{r syntax, eval=FALSE}
R_mixedregls(fo, indata, idname,  outdata,  outresults, save_def, 
    title1=NULL, title2=NULL, PNINT=NULL, RNINT=NULL, SNINT=NULL, 
    CONV=NULL, NQ=NULL,AQUAD=NULL, MAXIT=NULL, STD=NULL, MISS=NULL,
    NCOV=NULL) 
```

```{r example1}

# The data
indata <- read.table("nels23_data.txt", header=TRUE)

# A formula
fo <- formula(math ~ homew + sex + schtype | schtype | schtype  )

# Minimal input to run
R_mixregls(fo, indata, idname="school", 
     outdata="nels_example1.dat", 
     outresults="nels_example1.out", 
     save_def="nels_example1.def",
     NCOV=1)  
```
## Take A Look AT output

### Models fit and parameter estimates

```{r summary_of_models}
model_summaries<- readLines("nels_example1.out")
noquote(model_summaries)
```

### Just the Numbers

```{r just_numbers}
readLines("MIXREGLS.EST")
```

### Iteration history

```{r Iteration_history}
readLines("MIXREGLS.ITS")
```

### Standardized Residuals Model 1 & Random Effects

For the following, you can change the label on the columns and read them in as data frame or you can do something like what's below.  

```{r model1_std_residuals}
std1 <- readLines("MIXREGLS.RE1")
```

```{BLUPS and SEs}
std2 <- readLines("MIXREGLS.RE2")
```

# Example 2

In this example, I will set some of the options, one of which I recomend that is to give meaingful titles that will appead in the output.

```{r add_titles}
R_mixregls(fo, indata, idname="school", 
     outdata="nels_example1.dat", 
     outresults="nels_example1.out", 
     save_def="nels_example1.def",
     title1="NELS Example 1",
     title2="Illusration of adding titles")

# Table of output
model_summaries<- readLines("nels_example1.out")
noquote(model_summaries)
```

# A more Complex model

This is a random intercept version of Model 5 from lecture notes.

```{r center_ses_and_mean_ses}
schMses <- as.data.frame(aggregate(ses~school, data=indata, FUN="mean"))
names(schMses) <- c('school', 'schMses')
indata <- merge(indata,schMses, by=c('school'))
indata$schCses <- indata$ses - indata$schMses
```

This is similar to model 5 in lecture notes and Model_Building R example

```{r more_complex}
fo <- formula(math ~ homew + sex + schCses + schMses + schtype  | schtype  | schtype )
```

```{r fit_more_complex}
R_mixregls(fo, indata, idname="school", 
     outdata="nels_complex.dat", 
     outresults="nels_complex.out", 
     save_def="nels_complex.def",
     title1="NELS Example 1",
     title2="Illusration of adding titles")
```

```{r}
# Table of output
model_summaries<- readLines("nels_complex.out")
noquote(model_summaries)
```

# Setting BS to none
 A formula
 
```{r no_bs}
fo <- formula(math ~ homew + sex + schtype | none | schtype  )

R_mixregls(fo, indata, idname="school", 
     outdata="nels_no_bs.dat", 
     outresults="nels_no_bs.out", 
     save_def="nels_no_bs.def",
     title1="Model without BS predictors of Variance",
     title2="NELS data with 23 schools",
     NCOV=0)

# Table of output
model_no_bs<- readLines("nels_no_bs.out")
noquote(model_no_bs)
```

# Summary of Input to Function

Required input: 

* fo            formula with general form response ~ fixed | Between | Within
* indata        Your input data set
* idname        Name of id or group index
* outdat        Name of file with data that mixreglsb.exe will use 
* outresults    Name of file with mixregls summary
* save_def      Name of file to save model definition


Options and their defaults are 

* title1        The title in output (default is nothing). Recommended
* title2        The secondary title in outout (default is nothing). Recomemended 
* PNINT         =0 include intercept (default) for mean sub-model, =1 don't include 
* RNINT         =0 intercept (default) for beween sub-model, =1 don't include
* SNINT         =0 intercept (default) for within sub-model, =1 don't include  
* CONV          convergence criterin, default is 1E-6
* NQ            =number of quadrature points, default is 11
* AQUAD         =1 do adaptive quadrature, =0 don't
* MAXIT         maximum number of iterations, defaul is 200
* MISS          =0 no missing values ...see paper for how to include missing
* STD           =0 do not standardize covariates (default), =1 standardize
* NCOV          =1 linear association linear and random WS, for others see paper.
