7 MLM, Centering/Scaling: Reading Achievement (2-levels only)
library(tidyverse) # all things tidy
library(haven) # read in SPSS dataset
library(furniture) # nice table1() descriptives
library(stargazer) # display nice tables: summary & regression
library(texreg) # Convert Regression Output to LaTeX or HTML Tables
library(psych) # contains some useful functions, like headTail
library(performance) # ICC and R2 calculations
library(effects) # Effects for regression models
library(optimx) # Different optimizers to solve mlm's
library(lme4) # non-linear mixed-effects models
library(haven) # read in SPSS dataset
7.1 Background
The following example was included in the text “Multilevel Modeling in R” by Finch, Bolin, and Kelley (2016).
The datasets for this textbook may be downloaded from the website: http://www.mlminr.com/data-sets/.
I was unable to find any documentation on this dataset in the book or online, so I contacted the authors. There were unable to provide much either, but based on visual inspection designated the class of factor to thoes vairables that seem to represent categorical quantities. The labels for gender and class size are relative to the frequencies in the journal article the authors did point me to (although the samples sizes do not match up).
FOR THIS CHAPTER WE WILL IGNORE ALL LEVELS EXCEPT FOR STUDNETS BEING NESTED WITHIN SCHOOLS.
Read the SPSS data in with the haven
package .
<- haven::read_sav("http://www.mlminr.com/data-sets/Achieve.sav?attredirects=0") data_raw
Declare all categorical variables to be factors and apply labels where meaningful.
Student-specific
*gender
= Male or Female
*age
= Age, in months
*gevocab
= Vocabulary Score
*geread
= Reading Score
Class-specific
*classsize
= category of class’s size
School-specific
*senroll
= school enrollment
*ses
= school’s SES level
<- data_raw %>%
data_achieve ::mutate_at(vars(id, region, corp, school, class), factor) %>%
dplyr::mutate(gender = gender %>%
dplyrfactor(labels = c("Female", "Male"))) %>%
::mutate(classize = classize %>%
dplyrfactor(labels = c("12-17", "18-21",
"22-26", ">26"))) %>%
::select(id, region, corp, school, class, # Identifiers
dplyr# Pupil-level vars
gender, age, geread, gevocab, # Class-Level vars
classize, # School-level vars senroll, ses)
7.1.1 Sample Structure
It is obvious that the sample is hiarchical in nature. The nesting starts with students
(level 1) nested within class
(level 2), which are further nested within school
(level 3), corp
(level 4), and finally region
(level 5).
For this chapter we will only focus on TWO levels: students (level 1) are the units on which the outcome is measured and schools (level 2) are the units in which they are nested.
The number of regions = 9:
<- data_achieve %>%
num_regions ::group_by(region) %>%
dplyr::tally() %>%
dplyrnrow()
num_regions
[1] 9
The number of corps = 60:
<- data_achieve %>%
num_corps ::group_by(region, corp) %>%
dplyr::tally() %>%
dplyrnrow()
num_corps
[1] 60
The number of schools = 160
<- data_achieve %>%
num_schools ::group_by(region, corp, school) %>%
dplyr::tally() %>%
dplyrnrow()
num_schools
[1] 160
The number of classes = 568
<- data_achieve %>%
num_classes ::group_by(region, corp, school, class) %>%
dplyr::tally() %>%
dplyrnrow()
num_classes
[1] 568
The number of students = 10320
<- data_achieve %>% nrow
num_subjects
num_subjects
[1] 10320
7.2 Exploratory Data Analysis
7.2.1 Summarize Descriptive Statistics
<IP_ADDRESS> The stargazer
package
Most posters, journal articles, and reports start with a table of descriptive statistics. Since it tends to come first, this type of table is often refered to as Table 1. The stargazer()
function can be used to create such a table, but only for the entire dataset. I haven’t been able to find a way to get it to summarize subsamples and compare them in the standard format. Also, it only summarises continuous, not categorical variables.
# Knit to Website: type = "html"
# Knit to PDF: type = "latex"
# View on Screen: type = "text"
%>%
data_achieve ::select(classize, gender, geread, gevocab, age, ses, senroll) %>%
dplyrdata.frame() %>%
::stargazer(header = FALSE,
stargazertitle = "Summary of the numeric variables with `stargazer`",
type = "text")
| Statistic | N | Mean | St. Dev. | Min | Max |
| geread | 10,320 | 4.341 | 2.332 | 0.000 | 12.000 |
| gevocab | 10,320 | 4.494 | 2.368 | 0.000 | 11.200 |
| age | 10,320 | 107.529 | 5.060 | 82 | 135 |
| ses | 10,320 | 72.849 | 21.982 | 0.000 | 100.000 |
| senroll | 10,320 | 533.415 | 154.797 | 115 | 916 |
<IP_ADDRESS> The furniture
package
Tyson Barrett’s furniture package includes the extremely useful function table1()
which simplifies the common task of creating a stratified, comparative table of descriptive statistics. Full documentation can be accessed by executing ?furniture::table1
.
# Knit to Website: output = "html"
# Knit to PDF: output = "latex2"
# View on Screen: output = ""text", or "markdown", "html"
%>%
data_achieve :: table1("Reading Score" = geread,
furniture"Vocabulary Score" = gevocab,
"Age (in months)" = age,
"School SES" = ses,
"School's Enrollment" = senroll,
splitby = ~ gender, # var to divide sample by
test = TRUE, # test groups different?
caption = "Summary of the numeric variables with `table1`",
output = "html")
| Female | Male | P-Value | |
|---|---|---|---|
| n = 5143 | n = 5177 | ||
| Reading Score | 0.218 | ||
| 4.4 (2.3) | 4.3 (2.3) | ||
| Vocabulary Score | <.001 | ||
| 4.6 (2.4) | 4.4 (2.3) | ||
| Age (in months) | <.001 | ||
| 107.1 (5.0) | 107.9 (5.1) | ||
| School SES | 0.155 | ||
| 72.5 (22.3) | 73.2 (21.7) | ||
| School’s Enrollment | 0.483 | ||
| 532.3 (154.5) | 534.5 (155.1) |
7.2.2 Visualization of Raw Data
<IP_ADDRESS> Level One Plots: Disaggregate or ignore higher levels
For a first look, its useful to plot all the data points on a single scatterplot as displayed in the previous plot. Due to the large sample size, many points end up being plotted on top of or very near each other (overplotted). When this is the case, it can be useful to use geom_binhex()
rather than geom_point()
so the color saturation of the hexigons convey the number of points at that location, as seen in Figure \(\ref{fig:hexbin}\).
Note: I had to manually install the package hexbin
for the geom_hex()
to run.
%>%
data_achieve ggplot() +
aes(x = gevocab,
y = geread) +
stat_binhex(colour = "grey85", na.rm = TRUE) + # outlines
scale_fill_gradientn(colors = c("grey80","navyblue"), # fill color extremes
name = "Frequency", # legend title
na.value = NA) + # color for count = 0
theme_bw()
<IP_ADDRESS> Multilevel plots: illustrate two nested levels
Up to this point, all investigation of this dataset has been only at the pupil level and any nesting or clustering within schools has been ignored. Plotting is a good was to start to get an idea of the school-to-school variability. This figure displays four handpicked school to illustrate the degreen of school-to-school variability in the association between vocab and reading scores.
%>%
data_achieve ::filter(school %in% c(1321, 6181,
dplyr6197, 6823)) %>% # choose school numbers
ggplot(aes(x = gevocab,
y = geread))+
geom_count() + # creates points, size by overplotted number
geom_smooth(method = "lm") + # linear model (OLS)
facet_wrap(~ school) + # panels by school
theme_bw()
Another way to explore the school-to-school variability is to plot the linear model fit independently to each of the schools. This next figure displays only the smooth lines without the standard error bands or the raw data in the form of points or hexagons.
%>%
data_achieve ggplot(aes(x = gevocab,
y = geread)) +
geom_smooth(aes(group = school),
method = "lm",
se = FALSE, # do NOT want the SE bands
size = 0.3) +
geom_smooth(method = "lm",
se = FALSE,
color = "red", # do NOT want the SE bands
size = 2) + # make the lines thinner
theme_bw()
Due to the high number of schools, the figure with all the school’s independent linear regression lines resembles a hairball and is hard to deduce much about individual schools. By using the facet_grid()
layer, we can separate the schools out so better see school-to-school variability. It also allows investigation of higher level predictors, such as the school’s SES (median split with ntile(var, 2)
) and class size.
%>%
data_achieve ::mutate(ses2 = ntile(ses, 2) %>% # median split
dplyrfactor(labels = c("SES: Lower Half",
"SES: Upper Half"))) %>%
::mutate(senroll = ntile(senroll, 3) %>%
dplyrfactor(labels = c("Enroll: Smallest Third",
"Enroll: Middle Third",
"Enroll: Largest Third"))) %>%
ggplot(aes(x = gevocab,
y = geread,
group = school)) + # separates students into schools
geom_smooth(method = "lm",
se = FALSE,
size = 0.3,
color = "black",
alpha = .2) +
theme_bw() +
facet_grid(senroll ~ ses2) # makes separate panels (rows ~ columns)
7.3 Single-Level Regression
7.3.1 Fit Nested Models
Ignoring the fact that students are nested or clustered within schools, is called dissagregating. This treats all students as independent units.
# linear model - ignores school (for reference only)
<- lm(formula = geread ~ 1, # intercept only
fit_read_lm_0 data = data_achieve)
<- lm(formula = geread ~ gevocab , # one predictor
fit_read_lm_1 data = data_achieve)
<- lm(formula = geread ~ gevocab + age, # two predictors
fit_read_lm_2 data = data_achieve)
<- lm(formula = geread ~ gevocab*age, # interation+main effects
fit_read_lm_3 data = data_achieve)
Now compare the models:
::knitreg(list(fit_read_lm_0,
texreg
fit_read_lm_1,
fit_read_lm_2,
fit_read_lm_3),custom.model.names = c("Null",
"1 IV",
"2 IV",
"Interaction"),
caption = "OLS: Investigate Fixed, Pupil-level Predictors",
caption.above = TRUE,
single.row = TRUE)
| Null | 1 IV | 2 IV | Interaction | |
|---|---|---|---|---|
| (Intercept) | 4.34 (0.02)*** | 1.96 (0.04)*** | 3.19 (0.42)*** | 5.28 (0.87)*** |
| gevocab | 0.53 (0.01)*** | 0.53 (0.01)*** | 0.01 (0.19) | |
| age | -0.01 (0.00)** | -0.03 (0.01)*** | ||
| gevocab:age | 0.00 (0.00)** | |||
| R2 | 0.00 | 0.29 | 0.29 | 0.29 |
| Adj. R2 | 0.00 | 0.29 | 0.29 | 0.29 |
| Num. obs. | 10320 | 10320 | 10320 | 10320 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
Assess the significance of terms in the last ‘best’ model
summary(fit_read_lm_3)
Call:
lm(formula = geread ~ gevocab * age, data = data_achieve)
Residuals:
Min 1Q Median 3Q Max
-6.2069 -1.1250 -0.4362 0.6041 8.6476
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.282607 0.869769 6.074 1.3e-09 ***
gevocab 0.009154 0.189113 0.048 0.961394
age -0.030814 0.008066 -3.820 0.000134 ***
gevocab:age 0.004830 0.001759 2.746 0.006039 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.965 on 10316 degrees of freedom
Multiple R-squared: 0.2902, Adjusted R-squared: 0.29
F-statistic: 1406 on 3 and 10316 DF, p-value: < 2.2e-16
::r2(fit_read_lm_3) performance
# R2 for Linear Regression
R2: 0.290
adj. R2: 0.290
anova(fit_read_lm_3)
Analysis of Variance Table
Response: geread
Df Sum Sq Mean Sq F value Pr(>F)
gevocab 1 16224 16223.9 4202.2783 < 2.2e-16 ***
age 1 34 33.7 8.7356 0.003128 **
gevocab:age 1 29 29.1 7.5419 0.006039 **
Residuals 10316 39827 3.9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
7.3.2 Visualize the Interaction
::Effect(focal.predictors = c("gevocab", "age"), # chooses default values for
effectsmod = fit_read_lm_3) # continuous vars
gevocab*age effect
age
gevocab 82 95 110 120 140
0 2.755830 2.355243 1.893028 1.584884 0.9685969
3 3.971571 3.759370 3.514523 3.351291 3.0248284
6 5.187313 5.163497 5.136018 5.117699 5.0810599
8 5.997807 6.099582 6.217015 6.295304 6.4518809
10 6.808301 7.035667 7.298012 7.472909 7.8227019
::Effect(focal.predictors = c("gevocab", "age"), # chooses default values for
effectsmod = fit_read_lm_3) %>% # continuous vars
data.frame() %>%
mutate(age = factor(age)) %>% # must make a factor to separate lines
ggplot(aes(x = gevocab,
y = fit,
color = age)) +
geom_point() +
geom_line()
Here is a better version of the plot.
Age is in months, so we want multiples of 12 for good visualization
summary(data_achieve$age)/12 # divide by 12 to change months to years
Min. 1st Qu. Median Mean 3rd Qu. Max.
6.833 8.667 8.917 8.961 9.250 11.250
A good set set of illustrative ages could be: 7, 9, and 11:
c(7, 9, 11) * 12 # times by 12 to change years to months
[1] 84 108 132
::Effect(focal.predictors = c("gevocab", "age"),
effectsmod = fit_read_lm_3,
xlevels = list(age = c(84, 108, 132))) %>% # age is in months
data.frame() %>%
mutate(age_yr = factor(age/12)) %>% # it would be nice to plot age in years
ggplot(aes(x = gevocab,
y = fit,
color = age_yr,
linetype = age_yr)) +
geom_line(size = 1.25) +
theme_bw() +
labs(title = "Best Linear Model - Disaggregated Data (OLS)",
x = "Vocabulary Score",
y = "Reading Score",
linetype = "Age (yrs)",
color = "Age (yrs)") +
theme(legend.position = c(0.85, 0.2),
legend.key.width = unit(2, "cm"),
legend.background = element_rect(color = "black")) +
scale_linetype_manual(values = c("solid", "longdash", "dotted")) +
scale_x_continuous(breaks = seq(from = 0, to = 11, by = 2)) +
scale_y_continuous(breaks = seq(from = 0, to = 11, by = 1))
7.4 MLM - Step 1: Null Model, only fixed and random intercepts
A so called Empty Model only includes random intercepts. No independent variables are involved, other the grouping or clustering variable that designates how level 1 units are nested within level 2 units. For a cross-sectional study design this would be the grouping variables, where as for longitudinal or repeated measures designs this would be the subject identifier. This nested structure variable should be set to have class factor
.
7.4.1 Fit the Model
<- lme4::lmer(geread ~ 1 + (1|school),
fit_read_0ml data = data_achieve,
REML = FALSE) # fit via ML (not the default)
<- lme4::lmer(geread ~ 1 + (1|school) ,
fit_read_0re data = data_achieve,
REML = TRUE) # fit = REML (the default)
Compare the two models to OLS:
::knitreg(list(fit_read_lm_0,
texreg
fit_read_0ml,
fit_read_0re),custom.model.names = c("OLS",
"MLM-ML",
"MLM-REML"),
caption = "MLM: NULL Model,two estimation methods",
caption.above = TRUE,
single.row = TRUE)
| OLS | MLM-ML | MLM-REML | |
|---|---|---|---|
| (Intercept) | 4.34 (0.02)*** | 4.31 (0.05)*** | 4.31 (0.05)*** |
| R2 | 0.00 | ||
| Adj. R2 | 0.00 | ||
| Num. obs. | 10320 | 10320 | 10320 |
| AIC | 46270.34 | 46274.31 | |
| BIC | 46292.06 | 46296.03 | |
| Log Likelihood | -23132.17 | -23134.15 | |
| Num. groups: school | 160 | 160 | |
| Var: school (Intercept) | 0.39 | 0.39 | |
| Var: Residual | 5.05 | 5.05 | |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
Notice that the estimate for the intercept is nearly the same in the linear regression and intercept only models, but the standard errors are quite different. When there is clustering in sample, the result of ignoring it is under estimation of the standard errors and over stating the significance of associations. This table was made with the screenreg()
function in the self named package. I tend to prefer this display over stargazer()
.
7.4.2 Estimate the ICC
First, ask for the variance compenents:
::VarCorr(fit_read_0re) %>%
lme4print(comp = c("Variance", "Std.Dev"),
digits = 4)
Groups Name Variance Std.Dev.
school (Intercept) 0.3915 0.6257
Residual 5.0450 2.2461
::get_variance(fit_read_0re) insight
$var.fixed
[1] 0
$var.random
[1] 0.3915154
$var.residual
[1] 5.045008
$var.distribution
[1] 5.045008
$var.dispersion
[1] 0
$var.intercept
school
0.3915154
\[ \begin{align*} \text{schools} \rightarrow \; & \sigma^2_{u0} = 0.6257^2 = 0.392 \\ \text{students within schools} \rightarrow \; & \sigma^2_{e} = 2.2461^2 = 5.045 \\ \end{align*} \]
Intraclass Correlation (ICC) Formula \[ \overbrace{\rho}^{\text{ICC}} = \frac{\overbrace{\sigma^2_{u0}}^{\text{Random Intercept}\atop\text{Variance}}} {\underbrace{\sigma^2_{u0}+\sigma^2_{e}}_{\text{Total}\atop\text{Variance}}} \tag{Hox 2.9} \]
Then you can manually calculate the ICC.
0.392 / (0.392 + 5.045)
[1] 0.07209858
Or you can use the icc()
function in the performance
package.
::icc(fit_read_0re) performance
# Intraclass Correlation Coefficient
Adjusted ICC: 0.072
Unadjusted ICC: 0.072
Note: On page 45 (Finch, Bolin, and Kelley 2016), the authors substituted standard deviations into the formula, rather than variances. The mistake is listed on their webpage errata (http://www.mlminr.com/errata) and is repeated through the text.
7.5 MLM - Step 2: Add Lower-level explanatory variables, fixed, ML
Variance Component models (steps 2 and 3) - decompose the INTERCEPT variance into different variance compondents for each level. The regression intercepts are assumed to varry ACROSS the groups, while the slopes are assumed fixed (no random effects).
Fixed effects selection should come prior to random effects. You should use Maximum Likelihood (ML) estimation when fitting these models.
-
IF: only level 1 predictors and random intercepts are
incorporated
- Then: MLM \(\approx\) ANCOVA .
7.5.1 Add pupil’s vocab score as a fixed effects predictor
<- lme4::lmer(geread ~ gevocab + (1|school),
fit_read_1ml data = data_achieve,
REML = FALSE) # to compare fixed var sig
<- lme4::lmer(geread ~ gevocab + (1|school),
fit_read_1re data = data_achieve,
REML = TRUE) # for R-sq calcs
::knitreg(list(fit_read_0ml,
texreg
fit_read_1ml),custom.model.names = c("Null",
"w Pred"),
caption = "MLM: Investigate a Fixed Pupil-level Predictor",
caption.above = TRUE,
digits = 4)
| Null | w Pred | |
|---|---|---|
| (Intercept) | 4.3068*** | 2.0231*** |
| (0.0548) | (0.0492) | |
| gevocab | 0.5130*** | |
| (0.0084) | ||
| AIC | 46270.3388 | 43132.4318 |
| BIC | 46292.0643 | 43161.3991 |
| Log Likelihood | -23132.1694 | -21562.2159 |
| Num. obs. | 10320 | 10320 |
| Num. groups: school | 160 | 160 |
| Var: school (Intercept) | 0.3885 | 0.0987 |
| Var: Residual | 5.0450 | 3.7661 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
<IP_ADDRESS> Assess Significance of Effects
Likelihood Ratio Test (LRT)
Since models 0 and 1 are nested models, only differing by the the inclusion or exclusion of the fixed effects predictor gevocab
, AND both models were fit via Maximum Likelihood, we can compare the model fit may be compared via the Likilihood-Ratio Test (LRT). The Likelihood Ratio value (L. Ratio) is found by subtracting the two model’s -2 * logLik
or deviance
values. Significance is judged by the Chi Squared distribution, using the difference in the number of parameters fit as the degrees of freedom.
anova(fit_read_0ml, fit_read_1ml)
Data: data_achieve
Models:
fit_read_0ml: geread ~ 1 + (1 | school)
fit_read_1ml: geread ~ gevocab + (1 | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_0ml 3 46270 46292 -23132 46264
fit_read_1ml 4 43132 43161 -21562 43124 3139.9 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
What does the model look like?
::Effect(focal.predictors = c("gevocab"),
effectsmod = fit_read_1ml) %>%
data.frame() %>%
ggplot(aes(x = gevocab,
y = fit)) +
geom_ribbon(aes(ymin = lower,
ymax = upper),
alpha = .3) +
geom_line() +
theme_bw()
<IP_ADDRESS> Proportion of Variance Explained
Extract the variance-covariance estimates:
BL = BAseline: The Null Model (fit via REML)
::VarCorr(fit_read_0re) %>%
lme4print(comp = c("Variance", "Std.Dev"),
digits = 4)
Groups Name Variance Std.Dev.
school (Intercept) 0.3915 0.6257
Residual 5.0450 2.2461
\[ \sigma^2_{u0-BL} = 0.392 \\ \sigma^2_{e-BL} = 5.045 \]
MC = Model to Compare: Model with Predictor (fit via REML)
::VarCorr(fit_read_1re) %>%
lme4print(comp = c("Variance", "Std.Dev"),
digits = 4)
Groups Name Variance Std.Dev.
school (Intercept) 0.09978 0.3159
Residual 3.76647 1.9407
\[ \sigma^2_{u0-MC} = 0.100 \\ \sigma^2_{e-MC} = 3.766 \]
Level 1 \(R^2\) - Snijders and Bosker
Found on page 47 (Finch, Bolin, and Kelley 2016), the proportion of variance in the outcome explained by predictor on level one is given by:
Snijders and Bosker Formula - Level 1 \[ R^2_1 = 1 - \frac{\sigma^2_{e-MC} + \sigma^2_{u0-MC}} {\sigma^2_{e-BL} + \sigma^2_{u0-BL}} \]
Note: This formula also apprears in the Finch errata. The subscripts in the denominator of the fraction should be for model 0, not model 1. The formula is given correctly here. They did substitute in the correct values.
Calculate the value by hand:
1 - (0.100 + 3.766)/(0.392 + 5.045)
[1] 0.2889461
Or use the performance
package to help out: (Note it is using a difference method…)
::r2(fit_read_1re) performance
# R2 for Mixed Models
Conditional R2: 0.295
Marginal R2: 0.276
This means nearly 30% of the variance in reading scores, above and beyond that accounted for by school membership (i.e. school makeup or school-to-school variation), is attributable to vocabulary scores.
Level 1 \(R^2\) - Raudenbush and Bryk
Hox, Moerbeek, and Van de Schoot (2017) presents this formula on page 58 of chapter 2
Raudenbush and Bryk Approximate Formula - Level 1 \[ approx\; R^2_1 = \frac{\sigma^2_{e-BL} - \sigma^2_{e-MC}} {\sigma^2_{e-BL} } \tag{Hox 4.8} \]
Calculate the value by hand:
5.045 - 3.766) / 5.045 (
[1] 0.2535183
Although slightly different in value and meaning, this value also conveys that vocabulary scores are highly associated with reading scores.
Level 2 \(R^2\) - Snijders and Bosker Formula Extended
Snijders and Bosker Formula Extended - Level 2 \[ R^2_2 = 1 - \frac{\frac{\sigma^2_{e-MC}}{B} + \sigma^2_{u0-MC}} {\frac{\sigma^2_{e-BL}}{B} + \sigma^2_{u0-BL}} \]
\(B\) is the average size of the Level 2 units (schools). Technically, you should use the harmonic mean, but unless the clusters differ greatly in size, it doesn’t make a huge difference.
- Average sample cluster size
/ num_schools num_subjects
[1] 64.5
- Calculate by hand:
1 - ((3.766 / 64.5) + 0.100) / ((5.045 / 64.5) + 0.391)
[1] 0.6624428
This means that over two-thirds in school mean reading levels may be explained by their student’s vocabulary scores.
Level 2 \(R^2\) - Raudenbush and Bryk
Raudenbush and Bryk Approximate Formula - Level 2 \[ R^2_1 = \frac{\sigma^2_{u0-BL} - \sigma^2_{u0-MC}} {\sigma^2_{u0-BL} } \tag{Hox 4.9} \]
0.392 - 0.100)/(0.392) (
[1] 0.744898
Remeber that these ‘variance accounted for’ estimations are not as straight forwards as we would like.
7.5.2 Investigate More Level 1 Predictors
Part of investigating lower level explanatory variables, is checking for interactions between these variables. The interaction between fixed effects is also considered to be a fixed effect, so we need to employ Maximum Likelihood estimation to compare nested models.
<- lmer(geread ~ gevocab + age + (1 | school), # add main effect of age
fit_read_2ml data = data_achieve,
REML = FALSE)
<- lmer(geread ~ gevocab*age + (1 | school), # add interaction between vocab and age
fit_read_3ml data = data_achieve,
REML = FALSE)
::knitreg(list(fit_read_1ml,
texreg
fit_read_2ml,
fit_read_3ml),custom.model.names = c("Only Vocab",
"Both Main Effects",
"Interaction"),
caption = "MLM: Investigate Other Fixed Pupil-level Predictors",
caption.above = TRUE,
digits = 4)
| Only Vocab | Both Main Effects | Interaction | |
|---|---|---|---|
| (Intercept) | 2.0231*** | 3.0049*** | 5.1874*** |
| (0.0492) | (0.4172) | (0.8666) | |
| gevocab | 0.5130*** | 0.5121*** | -0.0279 |
| (0.0084) | (0.0084) | (0.1881) | |
| age | -0.0091* | -0.0294*** | |
| (0.0038) | (0.0080) | ||
| gevocab:age | 0.0050** | ||
| (0.0017) | |||
| AIC | 43132.4318 | 43128.8201 | 43122.5687 |
| BIC | 43161.3991 | 43165.0293 | 43166.0198 |
| Log Likelihood | -21562.2159 | -21559.4100 | -21555.2844 |
| Num. obs. | 10320 | 10320 | 10320 |
| Num. groups: school | 160 | 160 | 160 |
| Var: school (Intercept) | 0.0987 | 0.0973 | 0.0977 |
| Var: Residual | 3.7661 | 3.7646 | 3.7614 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
<IP_ADDRESS> Assess Significance of Effects
Likelihood Ratio Test (LRT)
anova(fit_read_1ml, fit_read_2ml, fit_read_3ml)
Data: data_achieve
Models:
fit_read_1ml: geread ~ gevocab + (1 | school)
fit_read_2ml: geread ~ gevocab + age + (1 | school)
fit_read_3ml: geread ~ gevocab * age + (1 | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_1ml 4 43132 43161 -21562 43124
fit_read_2ml 5 43129 43165 -21559 43119 5.6117 1 0.017841 *
fit_read_3ml 6 43123 43166 -21555 43111 8.2514 1 0.004072 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Not only is student’s age predictive of their reading level (I could have guessed that), but that age moderated the relationship between vocabulary and reading.
<IP_ADDRESS> Visulaize the Interation
Visulaizations are extremely helpful to interpred interactions.
summary(data_achieve$age)
Min. 1st Qu. Median Mean 3rd Qu. Max.
82.0 104.0 107.0 107.5 111.0 135.0
::Effect(focal.predictors = c("gevocab", "age"), # variables involved in the interaction
effectsmod = fit_read_3ml,
xlevels = list(age = c(84, 108, 132))) %>% # age is in months
data.frame() %>%
mutate(age_yr = factor(age/12)) %>% # it would be nice to plot age in years
ggplot(aes(x = gevocab,
y = fit,
color = age_yr)) +
geom_line() +
theme_bw()
There is a positive association between vocabulary and reading, but it is strongest for older childred. Among younger children, reading scores are more stable across vocabulary differences.
7.6 MLM - Step 3: Higher-level explanatory variables, fixed, ML
School enrollment (senroll
) applies to each school as a whole. When a variable is measured at a higher level, all units in the same group have the same value. In this case, all student in the same school have the same value for senroll
.
<- lme4::lmer(geread ~ gevocab*age + senroll + (1 | school),
fit_read_4ml data = data_achieve,
REML = FALSE)
::knitreg(list(fit_read_0ml,
texreg
fit_read_3ml,
fit_read_4ml),custom.model.names = c("Null",
"Level 1 only",
"Level 2 Pred"),
caption = "MLM: Investigate a Fixed School-Level Predictor",
caption.above = TRUE,
single.row = TRUE)
| Null | Level 1 only | Level 2 Pred | |
|---|---|---|---|
| (Intercept) | 4.31 (0.05)*** | 5.19 (0.87)*** | 5.24 (0.87)*** |
| gevocab | -0.03 (0.19) | -0.03 (0.19) | |
| age | -0.03 (0.01)*** | -0.03 (0.01)*** | |
| gevocab:age | 0.01 (0.00)** | 0.01 (0.00)** | |
| senroll | -0.00 (0.00) | ||
| AIC | 46270.34 | 43122.57 | 43124.31 |
| BIC | 46292.06 | 43166.02 | 43175.01 |
| Log Likelihood | -23132.17 | -21555.28 | -21555.16 |
| Num. obs. | 10320 | 10320 | 10320 |
| Num. groups: school | 160 | 160 | 160 |
| Var: school (Intercept) | 0.39 | 0.10 | 0.10 |
| Var: Residual | 5.05 | 3.76 | 3.76 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
7.6.1 Assess Significance of Effects
Likelihood Ratio Test (LRT)
anova(fit_read_0ml, fit_read_3ml, fit_read_4ml)
Data: data_achieve
Models:
fit_read_0ml: geread ~ 1 + (1 | school)
fit_read_3ml: geread ~ gevocab * age + (1 | school)
fit_read_4ml: geread ~ gevocab * age + senroll + (1 | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_0ml 3 46270 46292 -23132 46264
fit_read_3ml 6 43123 43166 -21555 43111 3153.7701 3 <2e-16 ***
fit_read_4ml 7 43124 43175 -21555 43110 0.2548 1 0.6137
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
School enrollment (or size) does not seem be related to reading scores.
7.7 MLM - Step 4: Explanatory variables predict Slopes, random, REML
Random Coefficient models - decompose the SLOPE variance BETWEEN groups.
The fixed effect of the predictor captures the overall association it has with the outcome (intercept), while the random effect of the predictor captures the group-to-group variation in the association (slope). Note: A variable can be fit as BOTH a fixed and random effect.
<- lme4::lmer(geread ~ gevocab*age + (1 | school), # refit the previous 'best' model via REML
fit_read_3re data = data_achieve,
REML = TRUE)
#fit_read_5re <- lmer(geread ~ gevocab + (gevocab | school),
# data = achieve,
# REML = TRUE) # failed to converge :(
<- lme4::lmer(geread ~ gevocab*age + (gevocab | school),
fit_read_5re data = data_achieve,
REML = TRUE,
control = lmerControl(optimizer = "optimx", # get it to converge
calc.derivs = FALSE,
optCtrl = list(method = "nlminb",
starttests = FALSE,
kkt = FALSE)))
::knitreg(list(fit_read_3re,
texreg
fit_read_5re),custom.model.names = c("Rand Int",
"Rand Int and Slopes"),
caption = "MLM: Investigate Random Effects",
caption.above = TRUE,
single.row = TRUE)
| Rand Int | Rand Int and Slopes | |
|---|---|---|
| (Intercept) | 5.19 (0.87)*** | 5.61 (0.87)*** |
| gevocab | -0.03 (0.19) | -0.14 (0.19) |
| age | -0.03 (0.01)*** | -0.03 (0.01)*** |
| gevocab:age | 0.01 (0.00)** | 0.01 (0.00)*** |
| AIC | 43155.49 | 43011.65 |
| BIC | 43198.95 | 43069.58 |
| Log Likelihood | -21571.75 | -21497.82 |
| Num. obs. | 10320 | 10320 |
| Num. groups: school | 160 | 160 |
| Var: school (Intercept) | 0.10 | 0.28 |
| Var: Residual | 3.76 | 3.66 |
| Var: school gevocab | 0.02 | |
| Cov: school (Intercept) gevocab | -0.06 | |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
<IP_ADDRESS> Assess Significance of Effect
Likelihood Ratio Test (LRT) for Random Effects
You can use the Chi-squared LRT test based on deviances even though
we fit our modesl with REML, since the models only differ in terms of
including/exclusing of a random effects; they have same fixed effects.
Just make sure to include the refit = FALSE
option.
anova(fit_read_3re, fit_read_5re, refit = FALSE)
Data: data_achieve
Models:
fit_read_3re: geread ~ gevocab * age + (1 | school)
fit_read_5re: geread ~ gevocab * age + (gevocab | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_3re 6 43155 43199 -21572 43143
fit_read_5re 8 43012 43070 -21498 42996 147.84 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There is evidence the effect child vocabulary has on reading varies across schools.
<IP_ADDRESS> Visualize the Model
What does the model look like?
::Effect(focal.predictors = c("gevocab", "age"),
effectsmod = fit_read_5re, # just different model
xlevels = list(age = c(84, 108, 132))) %>%
data.frame() %>%
::mutate(age_yr = factor(age/12)) %>%
dplyrggplot(aes(x = gevocab,
y = fit,
color = age_yr)) +
geom_line() +
theme_bw()
We are seeming much the same trends, but perhaps more separation between the lines.
7.8 MLM - Step 5: Cross-Level interactions between explanatory variables - fixed, ML
Cross-level interacitons involve variables at different levels. Here we will investigate the school-level enrollment moderating vocabulary’s effect since we say that vocab’s effect differs across schools (step 4).
Remember that an interaction beween fixed effects is also fixed.
<- lme4::lmer(geread ~ gevocab*age + (gevocab | school),
fit_read_5ml data = data_achieve,
REML = FALSE,
control = lmerControl(optimizer = "optimx",
calc.derivs = FALSE,
optCtrl = list(method = "nlminb",
starttests = FALSE,
kkt = FALSE)))
<- lme4::lmer(geread ~ gevocab*age + senroll + (gevocab | school),
fit_read_6ml data = data_achieve,
REML = FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
<- lme4::lmer(geread ~ gevocab*age + gevocab*senroll + (gevocab | school),
fit_read_7ml data = data_achieve,
REML = FALSE)
<- lme4::lmer(geread ~ gevocab*age*senroll + (gevocab | school),
fit_read_8ml data = data_achieve,
REML = FALSE)
If you get thelmer()
message: Some predictor
variables are on very different scales: consider rescaling, you
can trust your results, but you really should try re-scaling your
variables.
We are getting this message since gevoab
is on mostly a single digit scale,0 to 11.2, and age (in months) ranges in the low thripe-digits, 82 through 135, while school enrollment is in the mid-hundreds, 112-916. When we compute the interactions we get much, much larger values. Having variables on such widely different ranges of values can cause estimation problems.
%>%
data_achieve ::select(gevocab, age, senroll) %>%
dplyrsummary()
gevocab age senroll
Min. : 0.000 Min. : 82.0 Min. :115.0
1st Qu.: 2.900 1st Qu.:104.0 1st Qu.:438.0
Median : 3.800 Median :107.0 Median :519.0
Mean : 4.494 Mean :107.5 Mean :533.4
3rd Qu.: 5.200 3rd Qu.:111.0 3rd Qu.:644.0
Max. :11.200 Max. :135.0 Max. :916.0
For now, let us look at the results.
::knitreg(list(fit_read_5ml,
texreg
fit_read_6ml,
fit_read_7ml,
fit_read_8ml),custom.model.names = c("Level 1 only",
"Both Levels",
"Cross-Level",
"3-way"),
caption = "MLM: Investigate a Fixed Cross-level Interaction",
caption.above = TRUE,
digits = 3,
bold = TRUE)
| Level 1 only | Both Levels | Cross-Level | 3-way | |
|---|---|---|---|---|
| (Intercept) | 5.614*** | 5.603*** | 5.515*** | 5.503 |
| (0.872) | (0.878) | (0.894) | (3.111) | |
| gevocab | -0.138 | -0.138 | -0.114 | 0.530 |
| (0.192) | (0.192) | (0.198) | (0.661) | |
| age | -0.033*** | -0.033*** | -0.033*** | -0.033 |
| (0.008) | (0.008) | (0.008) | (0.029) | |
| gevocab:age | 0.006*** | 0.006*** | 0.006*** | 0.000 |
| (0.002) | (0.002) | (0.002) | (0.006) | |
| senroll | 0.000 | 0.000 | 0.000 | |
| (0.000) | (0.000) | (0.006) | ||
| gevocab:senroll | -0.000 | -0.001 | ||
| (0.000) | (0.001) | |||
| age:senroll | -0.000 | |||
| (0.000) | ||||
| gevocab:age:senroll | 0.000 | |||
| (0.000) | ||||
| AIC | 42979.711 | 42981.698 | 42983.452 | 42982.851 |
| BIC | 43037.646 | 43046.874 | 43055.870 | 43069.753 |
| Log Likelihood | -21481.855 | -21481.849 | -21481.726 | -21479.426 |
| Num. obs. | 10320 | 10320 | 10320 | 10320 |
| Num. groups: school | 160 | 160 | 160 | 160 |
| Var: school (Intercept) | 0.272 | 0.271 | 0.274 | 0.269 |
| Var: school gevocab | 0.019 | 0.019 | 0.019 | 0.019 |
| Cov: school (Intercept) gevocab | -0.062 | -0.062 | -0.063 | -0.061 |
| Var: Residual | 3.660 | 3.660 | 3.660 | 3.659 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
There is no evidence school enrollment moderates either of age or vocabulary’s effects.
<IP_ADDRESS> Assess Significance of Effects
Likelihood Ratio Test (LRT)
When you have a list of sequentially nested models, you can test them in order with one call to the anova()
funtion.
anova(fit_read_5ml, fit_read_6ml)
Data: data_achieve
Models:
fit_read_5ml: geread ~ gevocab * age + (gevocab | school)
fit_read_6ml: geread ~ gevocab * age + senroll + (gevocab | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_5ml 8 42980 43038 -21482 42964
fit_read_6ml 9 42982 43047 -21482 42964 0.0132 1 0.9087
anova(fit_read_5ml, fit_read_7ml)
Data: data_achieve
Models:
fit_read_5ml: geread ~ gevocab * age + (gevocab | school)
fit_read_7ml: geread ~ gevocab * age + gevocab * senroll + (gevocab | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_5ml 8 42980 43038 -21482 42964
fit_read_7ml 10 42983 43056 -21482 42963 0.2593 2 0.8784
anova(fit_read_5ml, fit_read_8ml)
Data: data_achieve
Models:
fit_read_5ml: geread ~ gevocab * age + (gevocab | school)
fit_read_8ml: geread ~ gevocab * age * senroll + (gevocab | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_5ml 8 42980 43038 -21482 42964
fit_read_8ml 12 42983 43070 -21479 42959 4.8597 4 0.302
7.9 Centering Predictors: Change Center
Centering variables measured on the lowest level only involves subtacting the mean from every value. The spread or standard deviation is not changed.
Although there are functions to automatically center and standardize variables, it is beneficial to manually create these variables, as it is more transparent and facilitates un-centering them later.
%>%
data_achieve ::select(gevocab, age, senroll) %>%
dplyrsummary()
gevocab age senroll
Min. : 0.000 Min. : 82.0 Min. :115.0
1st Qu.: 2.900 1st Qu.:104.0 1st Qu.:438.0
Median : 3.800 Median :107.0 Median :519.0
Mean : 4.494 Mean :107.5 Mean :533.4
3rd Qu.: 5.200 3rd Qu.:111.0 3rd Qu.:644.0
Max. :11.200 Max. :135.0 Max. :916.0
7.9.1 Use Centered Variables
NOTE: The models with CENTERED variables are able to be fit with the default optimizer settings and do not return the error: "unable to evaluate scaled gradientModel failed to converge: degenerate Hessian with 1 negative eigenvalues"
<- lme4::lmer(geread ~ I(gevocab - 4.494) * I(age - 107.5) +
fit_read_5ml_c I(gevocab - 4.494) | school),
(data = data_achieve,
REML = FALSE)
<IP_ADDRESS> Compare the Models
::knitreg(list(fit_read_5ml,
texreg
fit_read_5ml_c),custom.model.names = c("Raw Units",
"Grand Mean Centered"),
groups = list("Raw Scale" = 2:4,
"Mean Centered" = 5:7),
caption = "MLM: Investigate Centering Variables Involved in an Interaction",
caption.above = TRUE,
single.row = TRUE,
digits = 3)
| Raw Units | Grand Mean Centered | |
|---|---|---|
| (Intercept) | 5.614 (0.872)*** | 4.348 (0.033)*** |
| Raw Scale | ||
| gevocab | -0.138 (0.192) | |
| age | -0.033 (0.008)*** | |
| gevocab:age | 0.006 (0.002)*** | |
| Mean Centered | ||
| gevocab - 4.494 | 0.519 (0.014)*** | |
| age - 107.5 | -0.006 (0.004) | |
| I(gevocab - 4.494):age - 107.5 | 0.006 (0.002)*** | |
| AIC | 42979.711 | 42979.711 |
| BIC | 43037.646 | 43037.646 |
| Log Likelihood | -21481.855 | -21481.855 |
| Num. obs. | 10320 | 10320 |
| Num. groups: school | 160 | 160 |
| Var: school (Intercept) | 0.272 | 0.100 |
| Var: school gevocab | 0.019 | |
| Cov: school (Intercept) gevocab | -0.062 | |
| Var: Residual | 3.660 | 3.660 |
| Var: school I(gevocab - 4.494) | 0.019 | |
| Cov: school (Intercept) I(gevocab - 4.494) | 0.024 | |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
Notice that the interactions yield the exact same parameter estimates and significances, but the main effects (including the interactions) are different. Model fit statistics include \(-2LL\) are exactly the same, too.
::compare_performance(fit_read_5ml, fit_read_5ml_c) performance
# Comparison of Model Performance Indices
Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
--------------------------------------------------------------------------------------------------------------------------------
fit_read_5ml | lmerMod | 42979.7 (0.500) | 42979.7 (0.500) | 43037.6 (0.500) | 0.321 | 0.282 | 0.054 | 1.897 | 1.913
fit_read_5ml_c | lmerMod | 42979.7 (0.500) | 42979.7 (0.500) | 43037.6 (0.500) | 0.321 | 0.282 | 0.054 | 1.897 | 1.913
<IP_ADDRESS> Visualize the Model
What does the model look like?
First plot the model fit to the centered variables with all defaut settings.
::interact_plot(model = fit_read_5ml_c,
interactionspred = gevocab,
modx = age)
::interact_plot(model = fit_read_5ml_c,
interactionspred = gevocab,
modx = age,
modx.values = c(80, 110, 150),
interval = TRUE)
7.10 Rescaling Predictors: Change Units or Standardize
Where centering variables involved subtracting a set value, scalling a varaibles involves dividing by a set amount. When we both center to the mean and divide by the standard deviation, the new resulting varaible is said to be standardized (not to be confusing with normalizing, which is does not do).
To retain meaningful units, you can multiply or divide all the measured values of a variable by a set amount, like a multiple of 10. This retains the meaning behind the units while still bringing them into line with other variables in the model and can avoid some convergence issues.
7.10.1 Scale Varaibles
<IP_ADDRESS> Divide by a Meaningful Value
%>%
data_achieve ::select(gevocab, age, ses) %>%
dplyrsummary()
gevocab age ses
Min. : 0.000 Min. : 82.0 Min. : 0.00
1st Qu.: 2.900 1st Qu.:104.0 1st Qu.: 66.30
Median : 3.800 Median :107.0 Median : 81.70
Mean : 4.494 Mean :107.5 Mean : 72.85
3rd Qu.: 5.200 3rd Qu.:111.0 3rd Qu.: 87.80
Max. :11.200 Max. :135.0 Max. :100.00
For this situation, lets only divide SES by ten.
7.10.2 Use Scaled Variables
Using the new versions of our variables, investigate is SES has an effect, either in 2-way or 3-way interactions with age and vocabulary.
<- lme4::lmer(geread ~ I(gevocab - 4.494) * I(age - 107.5) +
fit_read_5ml_s I(gevocab - 4.494) | school),
(data = data_achieve,
REML = FALSE)
<- lme4::lmer(geread ~ I(gevocab - 4.494) * I(age - 107.5) +
fit_read_6ml_s I(gevocab - 4.494) * I(ses/10) +
I(gevocab - 4.494) | school),
(data = data_achieve,
REML = FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
<- lme4::lmer(geread ~ I(gevocab - 4.494) * I(age - 107.5) * I(ses/10) +
fit_read_7ml_s I(gevocab - 4.494) | school),
(data = data_achieve,
REML = FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
::knitreg(list(fit_read_5ml_s,
texreg
fit_read_6ml_s,
fit_read_7ml_s),custom.model.names = c("no SES",
"2-way",
"3-way"),
caption = "MLM: Investigate More Complex Fixed Interactions (with Scalling)",
caption.above = TRUE,
digits = 3)
| no SES | 2-way | 3-way | |
|---|---|---|---|
| (Intercept) | 4.348*** | 3.945*** | 3.942*** |
| (0.033) | (0.107) | (0.107) | |
| gevocab - 4.494 | 0.519*** | 0.671*** | 0.674*** |
| (0.014) | (0.051) | (0.051) | |
| age - 107.5 | -0.006 | -0.004 | -0.008 |
| (0.004) | (0.004) | (0.014) | |
| I(gevocab - 4.494):age - 107.5 | 0.006*** | 0.006** | 0.000 |
| (0.002) | (0.002) | (0.007) | |
| ses/10 | 0.057*** | 0.058*** | |
| (0.014) | (0.014) | ||
| I(gevocab - 4.494):ses/10 | -0.021** | -0.022** | |
| (0.007) | (0.007) | ||
| I(age - 107.5):ses/10 | 0.000 | ||
| (0.002) | |||
| I(gevocab - 4.494):I(age - 107.5):ses/10 | 0.001 | ||
| (0.001) | |||
| AIC | 42979.711 | 42928.813 | 42932.046 |
| BIC | 43037.646 | 43001.232 | 43018.948 |
| Log Likelihood | -21481.855 | -21454.407 | -21454.023 |
| Num. obs. | 10320 | 10320 | 10320 |
| Num. groups: school | 160 | 160 | 160 |
| Var: school (Intercept) | 0.100 | 0.079 | 0.079 |
| Var: school I(gevocab - 4.494) | 0.019 | 0.017 | 0.017 |
| Cov: school (Intercept) I(gevocab - 4.494) | 0.024 | 0.031 | 0.031 |
| Var: Residual | 3.660 | 3.662 | 3.662 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
anova(fit_read_5ml_s, fit_read_6ml_s, fit_read_7ml_s)
Data: data_achieve
Models:
fit_read_5ml_s: geread ~ I(gevocab - 4.494) * I(age - 107.5) + (I(gevocab - 4.494) | school)
fit_read_6ml_s: geread ~ I(gevocab - 4.494) * I(age - 107.5) + I(gevocab - 4.494) * I(ses/10) + (I(gevocab - 4.494) | school)
fit_read_7ml_s: geread ~ I(gevocab - 4.494) * I(age - 107.5) * I(ses/10) + (I(gevocab - 4.494) | school)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
fit_read_5ml_s 8 42980 43038 -21482 42964
fit_read_6ml_s 10 42929 43001 -21454 42909 54.8974 2 1.2e-12 ***
fit_read_7ml_s 12 42932 43019 -21454 42908 0.7674 2 0.6813
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There is evidence that SES moderates the main effect of vocabulary, after accounting for the interaction between age and vocabulary. But there is NOT evidence of a three-way interaction between vobaculary, age, and SES.
7.11 Final Model
Always refit the final model via REML.
<- lme4::lmer(geread ~ I(gevocab - 4.494) * I(age - 107.5) +
fit_read_6re_s I(gevocab - 4.494) * I(ses/10) +
I(gevocab - 4.494) | school),
(data = data_achieve,
REML = TRUE)
7.11.1 Table
::knitreg(fit_read_6re_s,
texregcaption = "MLM: Final Model",
caption.above = TRUE,
single.row = TRUE,
digits = 3)
| Model 1 | |
|---|---|
| (Intercept) | 3.945 (0.107)*** |
| gevocab - 4.494 | 0.671 (0.051)*** |
| age - 107.5 | -0.004 (0.004) |
| ses/10 | 0.057 (0.014)*** |
| I(gevocab - 4.494):age - 107.5 | 0.006 (0.002)** |
| I(gevocab - 4.494):ses/10 | -0.021 (0.007)** |
| AIC | 42976.511 |
| BIC | 43048.929 |
| Log Likelihood | -21478.255 |
| Num. obs. | 10320 |
| Num. groups: school | 160 |
| Var: school (Intercept) | 0.081 |
| Var: school I(gevocab - 4.494) | 0.018 |
| Cov: school (Intercept) I(gevocab - 4.494) | 0.032 |
| Var: Residual | 3.662 |
| ***p < 0.001; **p < 0.01; *p < 0.05 |
<IP_ADDRESS> Visualize the Model
Recall the scales that the revised variables are now on:
%>%
data_achieve ::select(gevocab, age, ses) %>%
dplyrsummary()
gevocab age ses
Min. : 0.000 Min. : 82.0 Min. : 0.00
1st Qu.: 2.900 1st Qu.:104.0 1st Qu.: 66.30
Median : 3.800 Median :107.0 Median : 81.70
Mean : 4.494 Mean :107.5 Mean : 72.85
3rd Qu.: 5.200 3rd Qu.:111.0 3rd Qu.: 87.80
Max. :11.200 Max. :135.0 Max. :100.00
::interact_plot(model = fit_read_6ml_s, # model name
interactionspred = gevocab, # x-axis, main independent variable (continuous, ordinal)
modx = age, # lines by moderator, another independent variable
mod2 = ses) # panel by 2nd moderator, another indep var
Define values for the moderator(s):
::interact_plot(model = fit_read_6ml_s,
interactionspred = gevocab,
modx = age,
mod2 = ses,
modx.values = c(90, 110, 130),
mod2.values = c(20, 55, 90),
interval = TRUE)
Swap the moderators:
::interact_plot(model = fit_read_6ml_s,
interactionspred = gevocab,
mod2 = age,
modx = ses,
mod2.values = c(90, 110, 130),
modx.values = c(20, 55, 90),
interval = TRUE)
7.11.2 Interpretation of model
There is strong evidence that higher vocabulary scores correlate with higher reading scores. This relationship is strongest in low SES schools and among older students. This relationship is also weaker in younger students and those attending high SES schools.
See: Nakagawa and Schielzeth (2013) for more regarding model \(R^2\)
From the performance::r2()
documentation, for mixed models:
- marginal r-squared considers only the variance of the fixed effects
- conditional r-squared takes both the fixed and random effects into account
::r2(fit_read_6re_s) performance
# R2 for Mixed Models
Conditional R2: 0.322
Marginal R2: 0.289
Over 32% of the variance in student reading is attributable vocab, ses, age, and school-to-school differences, R^2 = .322.
Helpful links:
http://maths-people.anu.edu.au/~johnm/r-book/xtras/mlm-ohp.pdf
http://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
http://web.stanford.edu/class/psych252/section_2015/Section_week9.html
https://www.r-bloggers.com/visualizing-generalized-linear-mixed-effects-models-part-2-rstats-lme4/


Could Haley’s role in Trump administration taint a 2024 bid?
by Associated PressCOLUMBIA, S.C. (AP) — With the launch of her new memoir, former U.N. Ambassador Nikki Haley has again ignited speculation she may be gearing up for a possible White House bid, perhaps as early as 2024.
But even though the former South Carolina governor left the Trump administration on her own terms, a rarity in a White House that’s seen its share of turmoil, some strategists say she could struggle to win back GOP traditionalists who aren’t fans of the polarizing president.
Haley has been navigating that balancing act since 2016, when there was initially no love lost between then-governor Haley and then-candidate Trump.
But in her two years at the United Nations, Haley managed to tread a careful path of occasionally speaking out against Trump while not directly drawing his ire.


Usage: GB / GB
ISP: ---- bendbroadband
If my usage goes over 99% 95% 90% 80% 1% , notify me every 48 hours 24 hours 4 hours 2 hours 10 minutes 5 minutes 1 minute .
Enable Debug? No Yes
Version: 0.7.1


The parse driver currently supports several clocks with different query mechanisms. In order for you to find a sample that might be similar to a clock you might want to integrate into parse I'll sum up the major features of the clocks (this information is distributed in the parse/clk_*.c and ntpd/refclock_parse.c files).
Meinberg: start=<STX>, end=<ETX>, sync on start pattern="\2D: . . ;T: ;U: . . ; \3" pattern="\2 . . ; ; : : ; \3" pattern="\2 . . ; ; : : ; : ; ; . . "
Meinberg is a German manufacturer of time code receivers. Those clocks have a pretty common output format in the stock version. In order to support NTP Meinberg was so kind to produce some special versions of the firmware for the use with NTP. So, if you are going to use a Meinberg clock please ask whether there is a special Uni Erlangen version. You can reach Meinberg via the Web. Information can also be ordered via eMail from<EMAIL_ADDRESS>General characteristics:
Meinberg clocks primarily output pulse per second and a describing ASCII string. This string can be produced in two modes: either upon the reception of a question mark or every second. NTP uses the latter mechanism. DCF77 AM clocks have a limited accuracy of a few milliseconds. The DCF77 PZF5xx variants provide higher accuracy and have a pretty good relationship between RS232 time code and the PPS signal. Except for early versions of the old GPS166 receiver type, Meinberg GPS receivers have a very good timing relationship between the datagram and the pulse. The beginning of the start bit of the first character has basically the same accuracy as the PPS signal, plus a jitter of up to 1 bit time depending on the selected baud rate, i.e. 52 μs @ 19200. PPS support should always be used, if possible, in order to yield the highest possible accuracy.
The preferred tty setting for Meinberg DCF77 receivers is 9600/7E2:
CFLAG (B9600|CS7|PARENB|CREAD|HUPCL) IFLAG (IGNBRK|IGNPAR|ISTRIP) OFLAG 0 LFLAG 0
The tty setting for Meinberg GPS16x/17x receivers is 19200/8N1:
CFLAG (B19200|CS8|PARENB|CREAD|HUPCL) IFLAG (IGNBRK|IGNPAR|ISTRIP) OFLAG 0 LFLAG 0
All clocks should be run at datagram once per second.
Format of the Meinberg standard time string:
<STX>D:dd.mm.yy;T:w;U:hh.mm.ss;uvxy<ETX> pos: 0<PHONE_NUMBER>11111111222222222233 3 1<PHONE_NUMBER>23456789012345678901 2 <STX> = start-of-text, ASCII code 0x02 dd.mm.yy = day of month, month, year of the century, separated by dots w = day of week (1..7, Monday = 1) hh:mm:ss = hour, minute, second, separated by dots u = '#' for GPS receivers: time is not synchronized '#' for older PZF5xx receivers: no correlation, not synchronized '#' for other devices: never sync'ed since powerup ' ' if nothing of the above applies v = '*' for GPS receivers: position has not been verified '*' for other devices: freewheeling based on internal quartz ' ' if nothing of the above applies x = 'U' if UTC time is transmitted 'S' if daylight saving time is active ' ' if nothing of the above applies y = '!' during the hour preceding start or end of daylight saving time 'A' during the hour preceding a leap second ' ' if nothing of the above applies <ETX> = end-of-text, ASCII code 0x03
Format of the Uni Erlangen time string for PZF5xx receivers:
<STX>dd.mm.yy; w; hh:mm:ss; tuvxyza<ETX> pos: 0<PHONE_NUMBER>11111111222222222233 3 1<PHONE_NUMBER>23456789012345678901 2 <STX> = start-of-text, ASCII code 0x02 dd.mm.yy = day of month, month, year of the century, separated by dots w = day of week (1..7, Monday = 1) hh:mm:ss = hour, minute, second, separated by colons t = 'U' if UTC time is transmitted, else ' ' u = '#' for older PZF5xx receivers: no correlation, not synchronized '#' for PZF511 and newer: never sync'ed since powerup ' ' if nothing of the above applies v = '*' if freewheeling based on internal quartz, else ' ' x = 'S' if daylight saving time is active, else ' ' y = '!' during the hour preceding start or end of daylight saving time, else ' ' z = 'A' during the hour preceding a leap second, else ' ' a = 'R' alternate antenna (reminiscent of PZF5xx), usually ' ' for GPS receivers <ETX> = end-of-text, ASCII code 0x03
Format of the Uni Erlangen time string for GPS16x/GPS17x receivers:
<STX>dd.mm.yy; w; hh:mm:ss; +uu:uu; uvxyzab; ll.lllln lll.lllle hhhhm<ETX> pos: 0<PHONE_NUMBER>111111112222222222333333333344444444445555555555666666 6 1<PHONE_NUMBER>234567890123456789012345678901234567890123456789012345 6 <STX> = start-of-text, ASCII code 0x02 dd.mm.yy = day of month, month, year of the century, separated by dots w = day of week (1..7, Monday = 1) hh:mm:ss = hour, minute, second, separated by colons +uu:uu = offset to UTC in hours and minutes, preceded by + or - u = '#' if time is not synchronized, else ' ' v = '*' if position has not been verified, else ' ' x = 'S' if daylight saving time is active, else ' ' y = '!' during the hour preceding start or end of daylight saving time, else ' ' z = 'A' during the hour preceding a leap second, else ' ' a = 'R' alternate antenna (reminiscent of PZF5xx), usually ' ' for GPS receivers b = 'L' during a leap second, i.e. if the seconds field is 60, else ' ' ll.lllln = position latitude in degrees, 'n' can actually be 'N' or 'S', i.e. North or South lll.lllle = position longitude in degrees, 'e' can actually be 'E' or 'W', i.e. East or West hhhh = position altitude in meters, always followed by 'm' <ETX> = end-of-text, ASCII code 0x03
Examples for Uni Erlangen strings from GPS receivers:
\x02 09.07.93; 5; 08:48:26; +00:00; ; 49.5736N 11.0280E 373m \x03 \x02 08.11.06; 3; 14:39:39; +00:00; ; 51.9828N 9.2258E 176m \x03
The Uni Erlangen formats should be used preferably. Newer Meinberg GPS receivers can be configured to transmit that format, for older devices there may be a special firmware version available.
For the Meinberg parse look into clk_meinberg.c
RAWDCF: end=TIMEOUT>1.5s, sync each char (any char),generate psuedo time codes, fixed format
direct DCF77 code input
In Europe it is relatively easy/cheap the receive the german time code transmitter DCF77. The simplest version to process its signal is to feed the 100/200ms pulse of the demodulated AM signal via a level converter to an RS232 port at 50Baud. parse/clk_rawdcf.c holds all necessary decoding logic for the time code which is transmitted each minute for one minute. A bit of the time code is sent once a second.
The preferred tty setting is: CFLAG (B50|CS8|CREAD|CLOCAL) IFLAG 0 OFLAG 0 LFLAG 0
From "Zur Zeit", Physikalisch-Technische Bundesanstalt (PTB), Braunschweig und Berlin, März 1989
Timecode transmission:
AM: time marks are send every second except for the second before the next minute mark time marks consist of a reduction of transmitter power to 25% of the nominal level the falling edge is the time indication (on time) time marks of a 100ms duration constitute a logical 0 time marks of a 200ms duration constitute a logical 1
see the spec. (basically a (non-)inverted psuedo random phase shift) encoding:
FM: Second Contents 0 - 10 AM: free, FM: 0 11 - 14 free 15 R - alternate antenna 16 A1 - expect zone change (1 hour before) 17 - 18 Z1,Z2 - time zone 0 0 illegal 0 1 MEZ (MET) 1 0 MESZ (MED, MET DST) 1 1 illegal 19 A2 - expect leap insertion/deletion (1 hour before) 20 S - start of time code (1) 21 - 24 M1 - BCD (lsb first) Minutes 25 - 27 M10 - BCD (lsb first) 10 Minutes 28 P1 - Minute Parity (even) 29 - 32 H1 - BCD (lsb first) Hours 33 - 34 H10 - BCD (lsb first) 10 Hours 35 P2 - Hour Parity (even) 36 - 39 D1 - BCD (lsb first) Days 40 - 41 D10 - BCD (lsb first) 10 Days 42 - 44 DW - BCD (lsb first) day of week (1: Monday -> 7: Sunday) 45 - 49 MO1 - BCD (lsb first) Month 50 MO10 - 10 Months 51 - 53 Y1 - BCD (lsb first) Years 54 - 57 Y10 - BCD (lsb first) 10 Years 58 P3 - Date Parity (even) 59 - usually missing (minute indication), except for leap insertion
Schmid clock: needs poll, binary input, end='\xFC', sync start
The Schmid clock is a DCF77 receiver that sends a binary time code at the reception of a flag byte. The contents if the flag byte determined the time code format. The binary time code is delimited by the byte 0xFC.
TTY setup is: CFLAG (B1200|CS8|CREAD|CLOCAL) IFLAG 0 OFLAG 0 LFLAG 0
The command to Schmid's DCF77 clock is a single byte; each bit allows the user to select some part of the time string, as follows (the output for the lsb is sent first).
Bit 0: time in MEZ, 4 bytes *binary, not BCD*; hh.mm.ss.tenths Bit 1: date 3 bytes *binary, not BCD: dd.mm.yy Bit 2: week day, 1 byte (unused here) Bit 3: time zone, 1 byte, 0=MET, 1=MEST. (unused here) Bit 4: clock status, 1 byte, 0=time invalid, 1=time from crystal backup, 3=time from DCF77 Bit 5: transmitter status, 1 byte, bit 0: backup antenna bit 1: time zone change within 1h bit 3,2: TZ 01=MEST, 10=MET bit 4: leap second will be added within one hour bits 5-7: Zero Bit 6: time in backup mode, units of 5 minutes (unused here)
Trimble SV6: needs poll, ascii timecode, start='>', end='<', query='>QTM<', eol='<'
Trimble SV6 is a GPS receiver with PPS output. It needs to be polled. It also need a special tty mode setup (EOL='<').
TTY setup is: CFLAG (B4800|CS8|CREAD) IFLAG (BRKINT|IGNPAR|ISTRIP|ICRNL|IXON) OFLAG (OPOST|ONLCR) LFLAG (ICANON|ECHOK)
Special flags are:
PARSE_F_PPSPPS - use CIOGETEV for PPS time stamping PARSE_F_PPSONSECOND - the time code is not related to the PPS pulse (so use the time code only for the second epoch) Timecode<PHONE_NUMBER>111111111122222222223333333 / char<PHONE_NUMBER>012345678901234567890123456 \ posn >RTMhhmmssdddDDMMYYYYoodnnvrrrrr;*xx< Actual ----33445566600112222BB7__-_____--99- Parse >RTM 1 ;* < Check
ELV DCF7000: end='\r', pattern=" - - - - - - - \r"
The ELV DCF7000 is a cheap DCF77 receiver sending each second a time code (though not very precise!) delimited by '`r'
Timecode YY-MM-DD-HH-MM-SS-FF\r FF&0x1 - DST FF&0x2 - DST switch warning FF&0x4 - unsynchronised
HOPF Funkuhr 6021 mit serieller Schnittstelle Created by F.Schnekenbuehl<EMAIL_ADDRESS>from clk_rcc8000.c Nortel DASA Network Systems GmbH, Department: ND250 A Joint venture of Daimler-Benz Aerospace and Nortel.
hopf Funkuhr 6021 used with 9600,8N1, UTC via serial line "Sekundenvorlauf" ON ETX zum Sekundenvorlauf ON dataformat 6021 output time and date transmit with control characters transmit evry second
Type 6021 Serial Output format
000000000011111111 / char<PHONE_NUMBER>01234567 \ position sABHHMMSSDDMMYYnre Actual C4110046231195 Parse s enr Check s = STX (0x02), e = ETX (0x03) n = NL (0x0A), r = CR (0x0D) A B - Status and weekday A - Status 8 4 2 1 x x x 0 - no announcement x x x 1 - Summertime - wintertime - summertime announcement x x 0 x - Wintertime x x 1 x - Summertime 0 0 x x - Time/Date invalid 0 1 x x - Internal clock used 1 0 x x - Radio clock 1 1 x x - Radio clock highprecision B - 8 4 2 1 0 x x x - MESZ/MEZ 1 x x x - UTC x 0 0 1 - Monday x 0 1 0 - Tuesday x 0 1 1 - Wednesday x 1 0 0 - Thursday x 1 0 1 - Friday x 1 1 0 - Saturday x 1 1 1 - Sunday
The Computime receiver sends a datagram in the following format every minute
Timestamp T:YY:MM:MD:WD:HH:MM:SSCRLF Pos<PHONE_NUMBER>012345678901 2 3<PHONE_NUMBER>111111111122 2 2 Parse T: : : : : : : \r\n T Startcharacter "T" specifies start of the timestamp YY Year MM Month 1-12 MD Day of the month WD Day of week HH Hour MM Minute SS Second CR Carriage return LF Linefeed
The WHARTON 400A Series clock is able to send date/time serial messages in 7 output formats. We use format 1 here because it is the shortest. We set up the clock to send a datagram every second. For use with this driver, the WHARTON 400A Series clock must be set-up as follows :
Programmable Selected Option No Option BST or CET display 3 9 or 11 No external controller 7 0 Serial Output Format 1 9 1 Baud rate 9600 bps 10 96 Bit length 8 bits 11 8 Parity even 12 E
WHARTON 400A Series output format 1 is as follows :
Timestamp STXssmmhhDDMMYYSETX Pos 0<PHONE_NUMBER>1234 0<PHONE_NUMBER>1111 STX start transmission (ASCII 0x02) ETX end transmission (ASCII 0x03) ss Second expressed in reversed decimal (units then tens) mm Minute expressed in reversed decimal hh Hour expressed in reversed decimal DD Day of month expressed in reversed decimal MM Month expressed in reversed decimal (January is 1) YY Year (without century) expressed in reversed decimal S Status byte : 0x30 + bit 0 0 = MSF source 1 = DCF source bit 1 0 = Winter time 1 = Summer time bit 2 0 = not synchronised 1 = synchronised bit 3 0 = no early warning 1 = early warning


Roc Nation: Jay-Z Hasn't Commented On Colin Kaepernick!
by illseedThere were rumors, but Roc Nation is refuting them all.
(AllHipHop Rumors) Over the weekend, Colin Kaepernick and the NFL had an old fashioned Mexican Stand-Off - with two entities! You know what that is, right?
Basically, nobody is giving an inch. Kap didn't give them what they demanded and he forced his demands to be met. At the end of the day, about six teams came to see the former football player practice.
But there was a side conversation...reports immediately came out to say that Jay-Z was "disappointed," "mad" or "pissed" at the man he once supported. But, was Jay-Z really upset? It does seem logical if the early rumors were right: that Jay-Z helped Kap get a tryout from the NFL. However, Roc Nation has shot down the reports that suggest Jay has say anything.
I am no expert, but it looks like Kap may still have his old skills in place! You think?
That's that! At the end of the day, this is what happened for fans. Thanks NFL and Kap!
At least Jay-Z isn't adding more gas to the fire! People are going crazy over this! I really wish they could have made it happen.


Version 2.3.8 (2012-04-12)
New Features
Changes
- Default value for T=50%
- Default value for V=100%
- Default value for C=1000
- Default value for Link Direction = Forward
- Default value for Weight-By = Degree
- Default value for Objective = Risk
Bugs fixed
Version 2.3.7 (2011-03-01)
New Features
Changes
Bugs fixed
- Betweenness will no longer be recalculated every time. It should now only be recalculated if there have been changes to the network structure.
Version 2.3.6 (2011-2-28)
New Features
Changes
- Changed the label in the NA Objective Graph Calculation UI area from "Unfixed" to "Variable".
- Threat and Vulnerability values can now be specified with 2 points of decimal precision.
Bugs fixed
- Fixed an issue that could cause the objective graph and the NA variable budget to be displaying different values.
Version 2.3.5 (2010-12-08)
New Features
Changes
- Changed bidirectional links error message for flow analysis to be more clear.
Bugs fixed
- When switching between Cascade and Flow simulations on the simulation tab on graphs that fail to meet the necessary requirements to perform flow analysis, displayed error messages will be cleared.
Version 2.3.4 (2010-11-15)
New Features
Changes
Bugs fixed
- Switched the X-Axis and Y-Axis text labels for the NA simulation graph.
Version 2.3.3 (2010-11-04)
New Features
Changes
- Changing the "Unfixed budget" on the calculation tab will now cause the Objective Graph to blank if the calculation button is not depressed.
- Changing the simulation type on the simulation tab will now cause the simulation graph to blank if the simulation button is not depressed.
Bugs fixed
- Fixed the allocation order for the stackleberg equations used in the NA calculation tab.
Version 2.3.2 (2010-10-26)
New Features
Changes
Bugs fixed
- Fixed an error that could cause XOR and AND gates to give incorrect results in FT calculations.
- Fixed an error with Betweenness calculations.
- Fixed default position of node edit window on smaller screens.
- Fixed an error caused by setting the elimination costs to zero in the FT calculations.
Version 2.3.1 (2010-10-21)
New Features
Changes
Bugs fixed
- When "none" is selected for a NA objective function, nodes will be ranked based on their current weight.
- Fixed a bug where opening a new file after working in another one would cause saves to overwrite the first file, rather than opening a new "Save As" window.
Version 2.3.0 (2010-10-20)
New Features
- Added Flow Consequence to analysis for simulation and calculation tabs.
- Added Consequence Reduction, Threat Maximization, and enhanced Risk Minimization to the
NA calculations.
- Added cascade failure simulation to the tools and graphs palettes.
- Monetary Units have been added to the MBRA Preferences menu. These add a text
label next to all monetary units for user reference. They are solely text-based, and
do not affect any calculations.
Changes
- Modified the graph so that at its default bar width, the bars will not label bars that have 100% height. This was done to prevent text label overlap due to the decreased width of the bars.
- Re-ordered the default order of NA data table columns.
- Calculations performed in the "Calculations" tab of the NA tools window will now
only update when the Calculate button in that window is on. This includes objective graph
updates and rankings on nodes in the main window.
- Contagiousness is now a "Weight By" option, rather than an "Objective Function".
Bugs fixed
- Fixed FT paths in data table.
Version 2.2.0 (2010-07-19)
New Features
- The Fault Tree data table now shows the path of each node.
Version 2.1.0 (2010-03-08)
New Features
- MBRA files are now saved compressed.
- Contagiousness objective function is now supported.
- Height weight function is now supported.
- NA Node links may now have direction.
- Original and reduced sum of the Objective Function is now displayed in the calculate pallet.
- New chart added that displays the effects of budget on the objective function.
Changes
Bugs fixed
- Calculated values are now filtered during load and save.
Version 2.0.4 (2009-12-11)
Bugs fixed
- Rank now sorts correctly in the property table.
Version 2.0.3 (2009-12-11)
Bugs fixed
- Reduced Vulnerability was incorrectly being set to zero with zero budget.
Version 2.0.2 (2009-12-10)
Bugs fixed
- Gamma value was being cached an not being set to zero when budget went to zero for
both Network Analysis and Fault Tree.
Version 2.0.1 (2009-12-09)
Changes
- Log reporter optimized for better performance.
Bugs fixed
- Workaround for a bug in QDoubleSpinBox, when trying to input a new value the
insertion point would get messed up resulting in a wrong number being inputted.
- Fixed stack overflow in Object Module caused by recursive delete.
- Placed guards on math.log function calls in Lua calculate scripts.
Version 2.0.0 (2009-10-30)
New Features
- OpenStreeMaps integration.
- Network Analysis updated with lots of new features.
- Added prevention budget to Network Analysis.
- Added betweenness calculations to Network Analysis.
- Added several objective functions to Network Analysis.
- Fault Tree calculations are new and improved.
- Added support for multiple Fault Trees.
- Added a way to attach a Fault Tree to a Network Analysis node.
- Added creating a Fault Tree from several Network Analysis nodes.
- Added copy and paste to Fault Tree.


| Home | Trees | Indices | Help |
|
|---|
|
|
The ``obvious`` module ====================== Use it to import very obvious functions. :Example: >>> from obvious import add >>> add(1, 1) 2 This is a subtitle ------------------- You can say so many things here ! You can say so many things here ! You can say so many things here ! You can say so many things here ! You can say so many things here ! You can say so many things here ! You can say so many things here ! You can say so many things here ! This is another subtitle ------------------------ Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.
|
|||
|
|
|||
__package__ = None
|
|
Adds two numbers and returns the result. This add two real numbers and return a real result. You will want to use this function in any place you would usually use the ``+`` operator but requires a functional equivalent. :param a: The first number to add :param b: The second number to add :type a: int :type b: int :return: The result of the addition :rtype: int :Example: >>> add(1, 1) 2 >>> add(2.1, 3.4) # all int compatible types work 5.5 .. seealso:: sub(), div(), mul() .. warnings:: This is a completly useless function. Use it only in a tutorial unless you want to look like a fool. .. note:: You may want to use a lambda function instead of this. .. todo:: Delete this function. Then masturbate with olive oil. |
| Home | Trees | Indices | Help |
|
|---|
| Generated by Epydoc 3.0.1 on Wed Mar 2 11:12:46 2016 | http://epydoc.sourceforge.net |


Settings
Error: You may only have 5 reminders.
Error: You may not have less than 1 reminder.
Reminder 1
Reminder 2
Reminder 3
Reminder 4
Reminder 5
Schedule
| Period | Class | Time |
|---|---|---|
| 1 | H Spanish II | 8:00-8:50 |
| 2 | H Web Dev | 8:55-9:45 |
| 3 | H Geometry | 9:50-10:40 |
| 4 | English | 10:45-11:35 |
| 5 | Lunch | 11:40-12:30 |
| 6 | H Economics | 12:35-1:25 |
| 7 | Driver's Ed | 1:30-2:20 |
| 8 | H Chemistry | 2:25-3:15 |


The importance of a dynamic, clear, and user-friendly website can mean the difference of either making that next important connection, or missing out on an oppertunity. Coming from a marketing and sales background, I am fully aware of how the rewards and repercussions of marketing and informational tools such as a website can have on a company as well as all types of organizations.
Possessing the highest understanding of current client-side web development languages and utilizing the latest tools allow me to not only fix and improve your current site, but build a user-friendly, functional, and dynamic website from the ground up. A website that you and your company can be proud of.
If you are interested in the services that I offer, or just want to say "Hi", feel free to send me a message.


Click a square to toggle it. Refresh page for a new card. Click the button below to save your card as a png.
 Find us on GitHub (make a PR or open an issue).
LIKES GORLS!!!
(FREE)
By Maia McCormick: http://bit.ly/gay_bingo


{{media.node.edge_media_preview_like.count | number: '1.'}}
likes
Quick, relevant and actionable instagram analytics.
*Instagram has closed several public endpoints in April 2018. Development in progress to account for this and crunch out detailed analytics.


Juice
Smoothie
Shakes
Ice Cream
Coffee&Cake
JUICES
Orange Bliss Rs 50
Exotica Mango Rs 30
Pineapple Maple Rs 30
Pear Rambo Rs 60
Carrot Juice Rs 70


Model PipeWithScalarField visualizes a pipe and a scalar field along the pipe axis. The latter is shown by mapping scalar field to color values with a color map and utilizing this color at the perimeter associated with the corresponding axis location. Typically the scalar field value is a temperature, but might be also another quantity. Predefined color maps are available from MultiBody.Visualizers.Colors.ColorMaps and can be selected via parameter "colorMap". A color map with the corresponding scalar field values can be exported as vector-graphics in svg-format with function MultiBody.Visualizers.Colors.colorMapToSvg. The position and orientation of the center of the circle at the left end of the pipe is defined via parameters "r_0" and "R", respectively. The pipe axis is oriented along the x-axis of the local coordinate system described by "R", see figure below:
The color coding is shown in the next figure. It was generated with MultiBody.Visualizers.Colors.colorMapToSvg using the following call:
colorMapToSvg(Modelica.Mechanics.MultiBody.Visualizers.Colors.ColorMaps.jet(), height=50, nScalars=6, T_max=100, caption="Temperature in C");
| Name | Description |
|---|---|
| colorMap |


grid Class
grid module - object for creating and manipulating grid widgets
Item Index
Methods
attachValidationListener
-
elem
Attaches an event listener to a specific grid cell for the 'validated' event. Overrides the normal 'blur' behavior for cell editing
Parameters:
-
elem
Object
buildAggregatesRow
-
gridData
-
gridId
Gets the aggregations - either passed in with metadata if server-side paging, or calculates the aggregations if client-side paging.
Returns:
buildHeaderAggregations
-
gridData
-
gridId
Builds the header's aggregations row if specified to be displayed at the top of the grid
copyGridWidth
-
gridElem
Copies the grid's column's widths to subsequent page data so that a consistent width is maintained.
Parameters:
-
gridElem
Object
create
-
gridData
-
gridElem
Exposed on the grid module. Called to create a grid widget.
createGridContent
-
gridData
-
gridElem
Creates the body of the grid. Sets the text in each cell and wires up the event listeners
createGridHeaders
-
gridData
-
gridElem
Creates the grid's headers given the column metadata passed to the create function
createGridInstanceMethods
-
gridElem
-
gridId
Adds publicly accessible methods to the DOM element that was passed to the 'create' function to initialize the grid widget
getInitialGridData
-
dataSource
-
callback
Abstraction for getting grid data so the 'create' function doesn't care whether it was passed in, or if there is a server call to get it.
initializeGrid
-
id
-
gridData
-
gridElem
Initializes an instance of the grid after retrieving the dataSource data. Sets the internal instance of the grid's data, calls to create the footer and content
Merge-Sort algorithm for grid data client-side sorting
-
data
-
sortObj
-
type
Parameters:
Returns:
saveCellEditData
-
input
On blur or successful validation if using the validator, removes the input from the grid cell, saves the data in the alteredData array and set a dirty flag on the grid dom element if the value changed
Parameters:
-
input
Object- The input element that was edited
saveCellSelectData
-
select
On blur or successful validation if using the validator, removes the input from the grid cell, saves the data in the alteredData array and set a dirty flag on the grid dom element if the value changed
Parameters:
-
select
Object- The select dom element


 AutoTS takes care of a few of these:
Pretransforming the data optimally for each model
Trying an assortment of different algorithms
Trying an assortment of hyperparameters for each algorithm
Underlying Process¶
AutoTS works in the following way at present:
It begins by taking long data and converting it to a wide dataframe with DateTimeIndex
An initial train/test split is generated where the test is the most recent data, of forecast_length
A random template of models is generated and tested on the initial train/test
Models consist of a pre-transformation step (fill na options, outlier removal options, etc), and algorithm (ie ETS) and model paramters (trend, damped, etc)
The top models (selected by a combination of metrics) are recombined with random mutations for n_generations
A handful of the best models from this process go to cross validation, where they are re-assessed on new train/test splits.
The best model in validation is selected as best_model and used in the
.predict()
method to generate forecasts.
A simple example¶
# also: _hourly, _daily, _weekly, or _yearly
from autots.datasets import load_monthly
df_long = load_monthly()
from autots import AutoTS
model = AutoTS(forecast_length=3, frequency='infer',
ensemble='simple', drop_data_older_than_periods=240,
max_generations=5, num_validations=2)
model = model.fit(df_long, date_col='datetime',
value_col='value', id_col='series_id')
# Print the name of the best model
print(model)
Import of data¶
There are two shapes/styles of pandas DataFrame
which are accepted.
The first is long data, like that out of an aggregated sales-transaction table containing three columns identified to .fit()
as date_col {pd.Datetime}, value_col {the numeric or categorical data of interest}, and id_col {id string, if multiple series are provided}
.
Alternatively, the data may be in a wide format where the index is a pandas.DatetimeIndex
, and each column is a distinct data series.
Some AutoTS
parameters are only applicable if long style data is provided, which is the recommended input shape.
If your data is already wide (one column for each value), to bring to a long format:
df_long = df_wide.melt(id_vars=['datetime_col_name'],
var_name='series_id', value_name='value')
You can tailor the process in a few ways…¶
The simplest way to improve accuracy is to increase the number of generations max_generations=15
. Each generation tries new models, taking additional time but improving the accuracy. The nature of genetic algorithms, however, means there is no consistent improvement for each generation, and large number of generations will often only result in minimal performance gains.
Another approach that may improve accuracy is to set ensemble='all'
. Ensemble parameter expects a single string, and can for example be 'simple,dist'
, or 'horizontal'
. As this means storing more details of every model, this takes more time and memory.
A handy parameter for when your data is expected to always be 0 or greater (such as unit sales) is to set no_negatives=True
. This forces forecasts to be greater than or equal to 0.
A similar function is constraint=2.0
. What this does is prevent the forecast from leaving historic bounds set by the training data. In this example, the forecasts would not be allowed to go above max(training data) + 2.0 * st.dev(training data)
, as well as the reverse on the minimum side. A constraint of 0
would constrain forecasts to historical mins and maxes.
Another convenience function is drop_most_recent=1
specifing the number of most recent periods to drop. This can be handy with monthly data, where often the most recent month is incomplete.
drop_data_older_than_periods
provides similar functionality but drops the oldest data to speed up the process on large datasets.
remove_leading_zeroes=True
is useful for data where leading zeroes represent a process which has not yet started.
When working with many time series, it can be helpful to take advantage of subset=100
. Subset specifies the interger number of time series to test models on, and can be useful with many related time series (1000’s of customer’s sales). Usually the best model on a 100 related time series is very close to that tested on many thousands (or more) of series. This speeds up the model process in these cases, but does not work with horizontal
ensemble types.
Subset takes advantage of weighting, more highly-weighted series are more likely to be selected. Weighting is used with multiple time series to tell the evaluator which series are most important. Series weights are assumed to all be equal to 1, values need only be passed in when a value other than 1 is desired. Note for weighting, larger weights = more important.
Validation and Cross Validation¶
Firstly, all models are initially validated on the most recent piece of data. This is done because the most recent data will generally most closely resemble the forecast future.
With very small datasets, there may be not be enough data for cross validation, in which case num_validations
may be set to 0. This can also speed up quick tests.
Cross validation helps assure that the optimal model is stable over the dynamics of a time series.
Cross validation can be tricky in time series data due to the necessity of preventing data leakage from future data points.
Here, two methods of cross validation are in place, 'even'
and ‘backwards'
.
Even cross validation slices the data into equal chunks. For example, num_validations=3
would split the data into equal, progressive thirds (less the original validation sample). The final validation results would then include four pieces, the results on the three cross validation samples as well as the original validation sample.
Backwards cross validation works backwards from the most recent data. First the most recent forecast_length samples are taken, then the next most recent forecast_length samples, and so on. This makes it more ideal for smaller or fast-changing datasets.
Seasonal validation is supplied as 'seasonal n'
ie 'seasonal 364'
. It trains on the most recent data as usual, then valdations are n
periods back from the datetime of the forecast would be. For example with daily data, forecasting for a month ahead, and n=364
, the first test might be on May 2020, with validation on June 2019 and June 2018, the final forecast then of June 2020.
Only a subset of models are based from initial validation to cross validation. The number of models is set such as models_to_validate=10
. If a float in 0 to 1 is provided, it is treated as a % of models to select. If you suspect your most recent data is not fairly representative of the whole, it would be a good idea to increase this parameter.
A more detailed example:¶
Here, we are forecasting the traffice along Interstate 94 between Minneapolis and St Paul in (lovely) Minnesota. This is a great dataset to demonstrate a recommended way of including external variables - by including them as time series with a lower weighting.
Here weather data is included - winter and road construction being the major influencers for traffic and will be forecast alongside the traffic volume. This carries information to models such as RollingRegression
, VARMAX
, and VECM
.
Also seen in use here is the model_list
.
from autots.datasets import load_hourly
df_long = load_hourly()
# all other series assumed to be weight of 1
weights_hourly = {'traffic_volume': 20}
model_list = ['ZeroesNaive', 'LastValueNaive', 'MedValueNaive', 'GLS',
'ETS', 'RollingRegression', 'UnobservedComponents', 'VECM']
from autots import AutoTS
model = AutoTS(forecast_length=73, frequency='infer',
prediction_interval=0.95, ensemble='simple',
max_generations=5, num_validations=2,
validation_method='even',
model_list=model_list, models_to_validate=15,
drop_most_recent=1)
model = model.fit(df_long, date_col='datetime',
value_col='value', id_col='series_id',
weights=weights_hourly)
prediction = model.predict()
forecasts_df = prediction.forecast
Probabilistic forecasts are available for all models, but in many cases are just data-based estimates in lieu of model estimates, so be careful.
upper_forecasts_df = prediction.upper_forecast
lower_forecasts_df = prediction.lower_forecast
Model Lists¶
By default, most available models are tried. For a more limited subset of models, a custom list can be passed in, or more simply, a string, one of 'probabilistic', 'multivariate', 'fast', 'superfast', or 'all'
.
As of 0.2.0
the following models are included:
['ZeroesNaive', 'LastValueNaive', 'AverageValueNaive', 'GLS',
'GLM', 'ETS', 'ARIMA', 'FBProphet', 'RollingRegression',
'GluonTS', 'SeasonalNaive', 'UnobservedComponents', 'VARMAX',
'VECM', 'DynamicFactor', 'TSFreshRegressor', 'MotifSimulation',
'WindowRegression', 'VAR', 'TensorflowSTS', 'TFPRegression',
'ComponentAnalysis']
On large multivariate series, TSFreshRegressor
, DynamicFactor
and VARMAX
can be impractically slow.
Deployment and Template Import/Export¶
Many models can be reverse engineered with (relative) simplicity outside of AutoTS by placing the choosen parameters into Statsmodels or other underlying package.
There are some advantages to deploying within AutoTS using a reduced starting template. Following the model training, the top models can be exported to a .csv
or .json
file, then on next run only those models will be tried.
This allows for improved fault tolerance (by relying not on one, but several possible models and underlying packages), and some flexibility in switching models as the time series evolve.
One thing to note is that, as AutoTS is still under development, template formats are likely to change and be incompatible with future package versions.
# after fitting an AutoTS model
example_filename = "example_export.csv" # .csv/.json
model.export_template(example_filename, models='best',
n=15, max_per_model_class=3)
# on new training
model = AutoTS(forecast_length=forecast_length,
frequency='infer', max_generations=0,
num_validations=0, verbose=0)
model = model.import_template(example_filename, method='only') # method='add on'
print("Overwrite template is: {}".format(str(model.initial_template)))
Metrics¶
There are a number of available metrics, all combined together into a ‘Score’ which evaluates the best model. The ‘Score’ that compares models can easily be adjusted by passing through custom metric weights dictionary. Higher weighting increases the importance of that metric, while 0 removes that metric from consideration. Weights should be 0 or positive numbers, and can be floats as well as integers. This weighting is not to be confused with series weighting, which effects how equally any one metric is applied to all the series.
metric_weighting = {'smape_weighting' : 10,
'mae_weighting' : 1,
'rmse_weighting' : 5,
'containment_weighting' : 1,
'runtime_weighting' : 0,
'spl_weighting': 1,
'contour_weighting': 0}
model = AutoTS(forecast_length=forecast_length,
frequency='infer',
metric_weighting=metric_weighting)
It is wise to usually use several metrics. I often find the best sMAPE model, for example, is only slightly better in sMAPE than the next place model, but that next place model has a much better MAE and RMSE.
Warning: weights are not fully balanced 1 - 1 - 1. As such it is usually best to place your favorite metric an order of magnitude or more above the others.
sMAPE
is generally the most versatile metric across multiple series, but doesn’t handle forecasts with lots of zeroes well.
SPL
is Scaled Pinball Loss and is the optimal metric for assessing upper/lower quantile forecast accuracies.
Containment
measures the percent of test data that falls between the upper and lower forecasts, and is more human readable than SPL.
Contour
is a unique measure. It is designed to help choose models which when plotted visually appear similar to the actual. As such, it measures the % of points where the forecast and actual both went in the same direction, either both up or both down, but not the magnitude of that difference. Does not work with forecast_length=1.
Installation and Dependency Versioning¶
pip install autots
Requirements:¶
Python >= 3.5
numpy
pandas
sklearn >= 0.20.0 (ColumnTransformer)
>= 0.23.0 (PoissonReg)
statsmodels
Of these, numpy and pandas are critical. Limited functionality should exist without scikit-learn. Full functionality should be maintained without statsmodels, albeit with fewer available models.
pip install autots['additional']
Optional Requirements¶
holidays
fbprophet
gluonts (requires mxnet)
mxnet (mxnet-mkl, mxnet-cu91, mxnet-cu101mkl, etc.)
tensorflow >= 2.0.0
lightgbm
xgboost
psutil
Experimental & Dev Requirements¶
tensorflow-probability
fredapi
tsfresh
Hardware Acceleration with Intel CPU and Nvidia GPU for Ubuntu/Windows¶
Download Anaconda or Miniconda.
(install Visual Studio if on Windows for C compilers)
If you have an Nvidia GPU, download NVIDIA CUDA and CuDNN.
Intel MKL is included with anaconda
and offers significant performance gain for Intel CPUs.
You can check if your system is using mkl with numpy.show_config()
.
If you are on an AMD CPU, you do not want to be using MKL, OpenBLAS is better.
On Linux ARM-based systems, apt-get/yum (rather than pip) installs of numpy/pandas may install faster compilations.
conda create -n timeseries python=3.8
conda activate timeseries
# for simplicity:
conda install anaconda
# elsewise:
conda install numpy scipy
conda install -c conda-forge scikit-learn
pip install statsmodels
# check the mxnet documentation for various flavors of mxnet available
pip install mxnet
pip install gluonts
conda update anaconda
conda install -c conda-forge fbprophet
pip install tensorflow
pip install tensorflow-probability
Caveats and Advice¶
Short Training History¶
How much data is ‘too little’ depends on the seasonality and volatility of the data.
Minimal training data most greatly impacts the ability to do proper cross validation. Set num_validations=0
in such cases.
Since ensembles are based on the test dataset, it would also be wise to set ensemble=None
if num_validations=0
.
Too much training data.¶
Too much data is already handled to some extent by 'context_slicer'
in the transformations, which tests using less training data.
That said, large datasets will be slower and more memory intensive, for high frequency data (say hourly) it can often be advisable to roll that up to a higher level (daily, hourly, etc.).
Rollup can be accomplished by specifying the frequency = your rollup frequency, and then setting the agg_func=np.sum
or ‘mean’ or other appropriate statistic.
Lots of NaN in data¶
Various NaN filling techniques are tested in the transformation. Rolling up data to a less-frequent frequency may also help deal with NaNs.
Adding regressors and other information¶
future_
regressor, to make it clear this is data that will be know with high certainy about the future.
Such data about the future is rare, one example might be number of stores that will be (planned to be) open each given day in the future when forecast sales.
Only a handful of models support adding regressors, and not all handle multiple regressors.
The recommended way to provide regressors is as a pd.Series/pd.Dataframe with a DatetimeIndex.
Don’t know the future? Don’t worry, the models can handle quite a lot of parallel time series, which is another way to add information. Additional regressors can be passed through as additional time series to forecast as part of df_long. Some models here can utilize the additional information they provide to help improve forecast quality. To prevent forecast accuracy for considering these additional series too heavily, input series weights that lower or remove their forecast accuracy from consideration.
Categorical Data¶
Categorical data is handled, but it is handled crudely. For example, optimization metrics do not currently include any categorical accuracy metrics. For categorical data that has a meaningful order (ie ‘low’, ‘medium’, ‘high’) it is best for the user to encode that data before passing it in, thus properly capturing the relative sequence (ie ‘low’=1, ‘medium’=2, ‘high’=3).
Custom and Unusual Frequencies¶
Data must be coercible to a regular frequency. It is recommended the frequency be specified as a datetime offset as per pandas documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#dateoffset-objects Some models will support a more limited range of frequencies.
Tested Frequencies¶
Models¶
Model |
Dependencies |
Optional Dependencies |
Probabilistic |
Multiprocessing |
GPU |
Experimental |
|---|---|---|---|---|---|---|
ZeroesNaive |
||||||
LastValueNaive |
||||||
AverageValueNaive |
True |
|||||
SeasonalNaive |
||||||
GLS |
statsmodels |
|||||
GLM |
statsmodels |
joblib |
||||
ETS |
statsmodels |
joblib |
||||
UnobservedComponents |
statsmodels |
|||||
ARIMA |
statsmodels |
True |
joblib |
|||
VARMAX |
statsmodels |
True |
||||
DynamicFactor |
statsmodels |
True |
||||
VECM |
statsmodels |
|||||
VAR |
statsmodels |
True |
||||
FBProphet |
fbprophet |
True |
joblib |
|||
GluonTS |
gluonts, mxnet |
True |
yes |
|||
RollingRegression |
sklearn |
lightgbm, tensorflow |
sklearn |
some |
||
WindowRegression |
sklearn |
lightgbm, tensorflow |
sklearn |
some |
||
MotifSimulation |
sklearn.metrics.pairwise |
True |
True |
|||
TensorflowSTS |
tensorflow_probability |
True |
yes |
True |
||
TFPRegression |
tensorflow_probability |
True |
yes |
True |
||
ComponentAnalysis |
sklearn |
True |
||||
TSFreshRegressor |
tsfresh, sklearn |
True |


*** Your Daily Free Stuff!
** In This Issue…
- Free Birds Eye Recipe Ready Item at Walmart
- FREE Medium Two-Topping Handmade Pan Pizza at Dominos During Baseball Season
- Free Sample of Schiff Move Free Ultra Omega
- Free Rachael Ray Nutrish for Cats Sample
- Free Jiffy Mix Recipe Book
- Free Sample of Truvia Natural Sweetener
- Search Sweet Free Stuff
** Free Birds Eye Recipe Ready Item at Walmart
Walmart has Birds Eye Recipe Ready items for $1. Use this $1 off Birds Eye Recipe Ready printable coupon to make it FREE!
The post Free Birds Eye Recipe Ready Item at Walmart appeared first on SweetFreeStuff.com.&
&
** FREE Medium Two-Topping Handmade Pan Pizza at Dominos During Baseball Season
*Update: A NO-HITTER HAS BEEN THROWN! Free pizza giveaway begins at September 2nd at 3pm ET! 20,000 registered MLB.com users are eligible for a free Dominos pizza beginning 3pm ET on Tuesday, 59/2/2014. Come back to MLB.com/dominos at that time to redeem. When a pitcher throws a no-hitter this season, visit MLB.com/Dominono and enter your [&
The post FREE Medium Two-Topping Handmade Pan Pizza at Dominos During Baseball Season appeared first on SweetFreeStuff.com.&
* Comments
- The first no-hitter couldn’t get through after filling in info. … by Night Owl
- WENT TO PAGE AND IT SAID OPPS YOU HAVE REACHED A WRONG PAGE by JOEL ZAMEN
- love that your guys are doing this =) by Rua Hasan
- Never able to find out how to sign up…tried but no one … by Debbie
- never could get it through. Let me put it information but kept … by Jane. Morse
- Plus 2 more…
&
** Free Sample of Schiff Move Free Ultra Omega
Move Free Ultra® Omega is powered by Omega 3s in 1 smaller pill instead of 2 large Glucosamine and Chondroitin tablets. &uest a free sample of Schiff Move Free Ultra Omega.
The post Free Sample of Schiff Move Free Ultra Omega appeared first on SweetFreeStuff.com.&
* Comments
- I’d love to try this product. Thank you. by Karen
- Would love one by cynthia
- I would like a free sample of this product. by amber zelko
- I would like to try this product, thanks by Thomas
&
** Free Rachael Ray Nutrish for Cats Sample
Free Rachael Ray Nutrish for Cats Sample
The post Free Rachael Ray Nutrish for Cats Sample appeared first on SweetFreeStuff.com.&
* Comments
- Would like a sample but can’t find the form by janice
- I would love a sample for my cat. by susan
- I would like samples…love the dog food by Leona Duffie
- Love this! Got my sample last week and my cat loved it!! by Amanda
- Me too Janet by Janice
- Plus 12 more…
&
** Free Jiffy Mix Recipe Book
The post Free Jiffy Mix Recipe Book appeared first on SweetFreeStuff.com.&
&
** Free Sample of Truvia Natural Sweetener
The post Free Sample of Truvia Natural Sweetener appeared first on SweetFreeStuff.com.&


Pvjs Highlighting Example - Basic
You can pass in a
highlights
property in your input arguments.
Pvjs Diagram Instance 1: WP525


WITNESS LIST State Cultural and Recreational Resources Committee February 25, 2003-2:00P HB 294 For: Billingsley, Randy (Self and Virginia Vancleave with Daughters Republic of Texas) Yates, Charles M. (Jan DeVault) Yates, Charles M. (Texian Legacy Assn.) Agency Overview On: Bryan, J.P. (Texas Historical Commission) Collen, Terry (Texas Historical Commission) Nau, John (Texas Historical Commission) Oaks, Larry (Texas Historical Commission) Rudd, Peggy D. (State Library & Archives Commission) Seidenberg, Edward (Texas State Library & Archives)


Skip to content
Posts
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
You can embed an R code chunk like this:
summary(cars)
## speed dist ## Min. : 4.0 Min. : 2.00 ## 1st Qu.:12.0 1st Qu.: 26.00 ## Median :15.0 Median : 36.00 ## Mean :15.4 Mean : 42.98 ## 3rd Qu.
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
You can embed an R code chunk like this:
summary(cars)
## speed dist ## Min. : 4.0 Min. : 2.00 ## 1st Qu.:12.0 1st Qu.: 26.00 ## Median :15.0 Median : 36.00 ## Mean :15.4 Mean : 42.98 ## 3rd Qu.