SMHS LinearModeling LMM - Revision history

Pineaumi: /* Footnotes */

2016-05-23T13:40:23Z

‎Footnotes

Pineaumi: /* Scientific Inference */

2016-05-23T13:39:20Z

‎Scientific Inference

Pineaumi: /* Scientific Inference */

2016-05-23T13:31:33Z

‎Scientific Inference

Pineaumi: /* Fixed and Random Effects */

2016-05-23T13:29:41Z

‎Fixed and Random Effects

Imoubara at 12:56, 23 May 2016

2016-05-23T12:56:21Z

Pineaumi: /* Fixed and Random Effects */

2016-05-20T17:58:56Z

‎Fixed and Random Effects

Pineaumi: /* Scientific Inference */

2016-05-20T17:56:26Z

‎Scientific Inference

Pineaumi: /* Scientific Inference */

2016-05-20T16:00:48Z

‎Scientific Inference

Pineaumi: /* Fixed and Random Effects */

2016-05-19T17:25:53Z

‎Fixed and Random Effects

Pineaumi: /* Fixed and Random Effects */

2016-05-19T17:21:16Z

‎Fixed and Random Effects

@@ Line 685: / Line 685: @@
 * <sup>12</sup> http://www.statmethods.net/input/importingdata.html
 * <sup>13</sup> http://wiki.socr.umich.edu/index.php/SOCR_Data_MLB_HeightsWeights
 ==Next See==

@@ Line 390: / Line 390: @@
 “… Age affects Weight ''X''<sup>2</sup>(1)=51.79, p< 10<sup>-12</sup>), increasing it (annually) by about 0.89 lb ± 0.1222 (standard error) …”
-We are using a Chi-Square test because Wilk’s Theorem   states that $-2log(\frac{Likelihood\_Null\_model(m0)}{Likelihood\_full\_model(m1)}$) approaches a Chi-Square distribution with degrees of freedom equal to the number of parameters that differ between the models (in this case, only “Age”, so ''df''=1).
+We are using a Chi-Square test because Wilk’s Theorem<sup>14</sup> states that $-2log(\frac{Likelihood\_Null\_model(m0)}{Likelihood\_full\_model(m1)}$) approaches a Chi-Square distribution with degrees of freedom equal to the number of parameters that differ between the models (in this case, only “Age”, so ''df''=1).
 This likelihood-based inference approach is somewhat different from the classical t-tests, ANOVA and linear model inference. Instead of obtaining a p-value directly from the MLE model, we need to run a second-order analysis (in this case ANOVA on the 2 models) to compare them and derive a p-value.

@@ Line 676: / Line 676: @@
 <b>lmer.model.2 = lmer(Weight ~ Height + Age + <u>(1+Height|Team) + (1+Height|Position)</u>, data=data, REML=FALSE)</b>
-<sup>9</sup> http://www.rstudio.com/products/rstudio/download/
+===Footnotes===
-<sup>10</sup> http://cran.r-project.org/web/packages/lme4/lme4.pdf
+* <sup>7</sup> http://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_020108_HeightsWeights
+* <sup>8</sup> http://wiki.socr.umich.edu/index.php/SOCR_Data_MLB_HeightsWeights
-<sup>11</sup> http://www.r-tutor.com/r-introduction/data-frame/data-import
+* <sup>9</sup> http://www.rstudio.com/products/rstudio/download/
+* <sup>10</sup> http://cran.r-project.org/web/packages/lme4/lme4.pdf
-<sup>12</sup> http://www.statmethods.net/input/importingdata.html
+* <sup>11</sup> http://www.r-tutor.com/r-introduction/data-frame/data-import
+* <sup>12</sup> http://www.statmethods.net/input/importingdata.html
-<sup>13</sup> http://wiki.socr.umich.edu/index.php/SOCR_Data_MLB_HeightsWeights
+* <sup>13</sup> http://wiki.socr.umich.edu/index.php/SOCR_Data_MLB_HeightsWeights
 ==Next See==

@@ Line 9: / Line 9: @@
 ===Fixed and Random Effects===
-Linear models express relationships between data elements (variables) in terms of a (linear) function. For example, we can model weight as a function of height.
+Linear models express relationships between data elements (variables) in terms of a (linear) function. For example, we can model weight as a function of height<sup>7</sup>.
 Weight ~ Height + ε
@@ Line 15: / Line 15: @@
 Here, “height” is a fixed effect, and ε is our “error term” representing the deviations between the model predictions and observations (of weight) due to “random” factors that we cannot control experimentally. This tem (ε) is the “probabilistic” or “stochastic” part of the model. Let’s try to unpack “ε” and add complexity to it. In mixed (fixed and random effect) models, everything in the “systematic” part of your model works just like with linear models. If we change the random aspect of our model this leaves the systematic part (height) unchanged.
-Suppose we’re looking at the Baseball data   and try to identify a relationship that looks like this:
+Suppose we’re looking at the Baseball data<sup>8</sup> and try to identify a relationship that looks like this:
 Weight ~ position + ε

@@ Line 4: / Line 4: @@
 <B>Questions:</B>
-<li>What happens if data are not independent and identically distributed (IIDs)?</li>
+*What happens if data are not independent and identically distributed (IIDs)?
-<li>How to model multiple observations for the case/subject (across time, conditions, etc.)?</li>
+*How to model multiple observations for the case/subject (across time, conditions, etc.)?
 ===Fixed and Random Effects===

@@ Line 112: / Line 112: @@
 To add an additional random effect for position we expand the model:
-<blockquote>Weight ~ Height + position + (1|player) + (1|position) + ε</blockquote>
+Weight ~ Height + position + (1|player) + (1|position) + ε
 Note that a similar argument may be made for “Team”, there could be error-term dependencies or patterns reflecting the between-team random effects on weight – some teams may be bulkier than others. This is really well documented in sports (e.g., top notch Serbian and Spanish water polo teams have a vastly different physiques).
@@ Line 251: / Line 250: @@
 Adding Age as an additional fixed effect:
-<blockquote><b>lmer.model.2 = lmer(Weight ~ Height + Age + (1|Team) + (1|Position), data=data)</b></blockquote>
+<b>lmer.model.2 = lmer(Weight ~ Height + Age + (1|Team) + (1|Position), data=data)</b>
 This expands the original model object lmer.model with an age predictor. “Age” is added as a fixed effect, because the relationship between Age and Weight is systematic and predictable (not random). Older players are expected to have higher weights.

@@ Line 428: / Line 428: @@
 Let’s inspect the coefficients of the mixed model by Team and by Position:
-coef<b>(lmer.model.1)</b>
+ coef<b>(lmer.model.1)</b>
 <center>
@@ Line 530: / Line 530: @@
 As expected from a mixed-effect model, each Team and each Position is assigned a different intercept, as we specified the model with “(1|Team)” and “(1|Position)” to take into account by-Team and by-Position variability:
-<blockquote><b>lmer.model.1 <- lmer(Weight ~ Height + Age + (1|Team) + (1|Position), data=data, REML=FALSE)</b></blockquote>
+<b>lmer.model.1 <- lmer(Weight ~ Height + Age + (1|Team) + (1|Position), data=data, REML=FALSE)</b>
 The fixed effects <b>(Height + Age)</b> are all the same for all Team and Positions, as this model is a random intercept model where we account for baseline-differences in Weight assuming whatever the effect of Height is, it’s the same for all Teams and Positions.
@@ Line 544: / Line 544: @@
 Examining the coefficients of this updated model:
-coef<b>(lmer.model.2)</b>
+ coef<b>(lmer.model.2)</b>
 Generates the following output:
@@ Line 652: / Line 652: @@
 Next, we need to address the significance of these effects (i.e., compute the corresponding p-values using the likelihood ration test). We keep our model from above <b>(lmer.model.2)</b> and compare it to a new null model <b>(lmer.model.0)</b> using the likelihood ratio test:
-<blockquote><b>lmer.model.0 = lmer(Weight ~  Age + <U>(1+Height|Team) + (1+Height|Position)</U>, data=data, REML=FALSE)</b></blockquote>
+<b>lmer.model.0 = lmer(Weight ~  Age + <U>(1+Height|Team) + (1+Height|Position)</U>, data=data, REML=FALSE)</b>
-<blockquote><b>lmer.model.2 = lmer(Weight ~ <span style="background-color: #FFFF00">Height</span> + Age + <U>(1+Height|Team) + (1+Height|Position)</U>, data=data, REML=FALSE)</b></blockquote>
+<b>lmer.model.2 = lmer(Weight ~ <span style="background-color: #FFFF00">Height</span> + Age + <U>(1+Height|Team) + (1+Height|Position)</U>, data=data, REML=FALSE)</b>
 The null model <b>(lmer.model.0)</b> needs to have the same random effects structure as the full model <b>(lmer.model.2).</b> That is, all random slope parameters included in the full model must be present in the null model.
@@ Line 660: / Line 660: @@
 Let’s now do the likelihood ratio test:
-anova<b>(lmer.model.0, lmer.model.2)</b>
+ anova<b>(lmer.model.0, lmer.model.2)</b>
   Models:
@@ Line 675: / Line 675: @@
 <span style="background-color: #FFFF00">The full model <b>(lmer.model.2)</b> is focused on baseball players’ Weight. As Age is not a variable we are interested in studying, albeit it may affect Weight, we need to control for Age. Thus <b>lmer.model.</b>2 has random slopes for the effect of <b><u>Height</u></b> (by Team and Position) but not the effect of Age. In other words, we only modeled by-Team and by-Position variability in how Height affects Weight.</span>
-<blockquote><b>lmer.model.2 = lmer(Weight ~ Height + Age + <u>(1+Height|Team) + (1+Height|Position)</u>, data=data, REML=FALSE)</b></blockquote>
+<b>lmer.model.2 = lmer(Weight ~ Height + Age + <u>(1+Height|Team) + (1+Height|Position)</u>, data=data, REML=FALSE)</b>
 <sup>9</sup> http://www.rstudio.com/products/rstudio/download/

@@ Line 399: / Line 399: @@
 What happens if we compare the following two models:
-<blockquote>m0: Weight ~ 1		<span style="padding-left:40px"># intercept only constant model estimating the mean only</span></blockquote>
+m0: Weight ~ 1                  # intercept only constant model estimating the mean only<BR>
-<blockquote>m1: Weight ~ Height + Age	<span style="padding-left:20px"># Intercept and 2 covariates model</span></blockquote>
+m1: Weight ~ Height + Age	# Intercept and 2 covariates model
 Suppose the likelihood-ration test suggest significant differences between the 2 models, this would suggest that “m1” (full) and “m0” (null) are significantly different from one another. However, it would not attribute this difference to “Height” or “Age” alone. We may not be able to tease out conclusively which one of the two variables (or both) was crucial.
@@ Line 406: / Line 406: @@
 How about the possibility of having an Age-by-Height interaction effect? For instance, “Height” effect on Weight may be modulated through “Age”. If such inter-dependence between two factors (an interaction) is present or suspected, its impact may be assessed the same way:
-<blockquote>full model:		<span style="padding-left:20px">Weight ~ Height * Age</span></blockquote>
+full model: Weight ~ Height * Age<BR>
-<blockquote>reduced model:	<span style="padding-left:20px">Weight ~ Height + Age</span></blockquote>
+reduced model: Weight ~ Height + Age
 Interactions between two factors are commonly specified with a “*” rather than a “+”. Comparing these models via the likelihood ratio test and the anova() function yields a p-value quantifying the significance of the interaction term (Height * Age). A significant result implies that Height and Age are strongly inter-dependent on each other. An insignificant result implies that there is no strong evidence in these data for inter-dependence.
-<b>Note:</b>
+<b>Note:</b><BR>
-<blockquote>lm.1 <- lm(y ~ x:z) 		<span style="padding-left:20px"># includes only the interaction (colon, “:”) between x & z,</span></blockquote>
+lm.1 <- lm(y ~ x:z) 	# includes only the interaction (colon, “:”) between x & z,
-<blockquote><span style="padding-left:20px"># but not the original variables x and z</span></blockquote>
+&#35; but not the original variables x and z <BR>
-<blockquote>lm.2 <- lm(y ~ x + z + x:z) 	<span style="padding-left:20px"># this is a complete model for y using x & z</span></blockquote>
+lm.2 <- lm(y ~ x + z + x:z) # this is a complete model for y using x & z
-<blockquote>#In regressions with a lot of variables, the product (*) notation provides a shortcut to minimize typing.</blockquote>
+&#35;In regressions with a lot of variables, the product (*) notation provides a shortcut to minimize typing.<BR>
-<blockquote># We can indicate this using a simple multiplication symbol (*):</blockquote>
+&#35; We can indicate this using a simple multiplication symbol (*):<BR>
-<blockquote>lm.3 <- lm(y ~ x*z)		<span style="padding-left:20px"># equivalent to model lm.2, but shorter.</span></blockquote>
+lm.3 <- lm(y ~ x*z)	# equivalent to model lm.2, but shorter.
 Using this code, R will automatically include both variables in the regression in addition to the interaction between the two
@@ Line 428: / Line 428: @@
 Let’s inspect the coefficients of the mixed model by Team and by Position:
-<blockquote>coef<b>(lmer.model.1)</b></blockquote>
+coef<b>(lmer.model.1)</b>
 <center>
@@ Line 538: / Line 538: @@
 Thus, we may need a random slope model, where Team and Positions may have differing intercepts as well as different slopes for the Weight. In R this can be coded as:
-<blockquote><b>lmer.model.2 = lmer(Weight ~ Height + Age + (1+Height|Team) + (1+Height|Position), data=data, REML=FALSE)</b></blockquote>
+<b>lmer.model.2 = lmer(Weight ~ Height + Age + (1+Height|Team) + (1+Height|Position), data=data, REML=FALSE)</b>
 This new model only includes more complicated random effects. In this case, the notation “(1+Height|Team)” means that the model expects differing baseline-levels of Weight (the intercept, represented by 1) <b>and</b> differing levels of “Height”. The same is true for Position.
@@ Line 544: / Line 544: @@
 Examining the coefficients of this updated model:
-<blockquote>coef<b>(lmer.model.2)</b></blockquote>
+coef<b>(lmer.model.2)</b>
 Generates the following output:

@@ Line 213: / Line 213: @@
 The ''standard deviation column'' represents a measure of the variability for each random effect included in the model. “Team” has much less variability than Position. The <u>Residual row represents the variability that’s not due to either Team or Position.</u> This is the ε term – random deviations from the predicted values that are not due to Team or Position. This reflects the fact that body-size has some factors that affect Weight that are outside of the observations we have in this experiment. Similarly, to the fixed effect models, the fixed effects output in this mixed model output includes estimates, SE and T-statistic.
-The coefficient of “Height” (4.6973) is the slope for this variable. Finally, the Correlation of the Fixed Effects (H vs. Intercept) is reported to be negative! <u>The output under "''correlation of fixed effects''" has a different interpretation from the intuitive meaning.</u> It is not about the correlation of the variables, but about the <u>“expected correlation of the regression coefficients”.</u> This may, or may not, be related to multicollinearity. In this case, ''corr(H,Intrc)''= -0.991 suggests that <u>if we redo the experiment again a decrease of the Height coefficient will increase the intercept</u>, and vice-versa, an increase of the Height coefficient is expected to drive down the intercept term (which is currently =-142.8215). <mark><-check for proofreading, clarity</mark>
+The coefficient of “Height” (4.6973) is the slope for this variable. Finally, the Correlation of the Fixed Effects (H vs. Intercept) is reported to be negative! <u>The output under "''correlation of fixed effects''" has a different interpretation from the intuitive meaning.</u> It is not about the correlation of the variables, but about the <u>“expected correlation of the regression coefficients”.</u> This may, or may not, be related to multicollinearity. In this case, ''corr(H,Intrc)''= -0.991 suggests that <u>if we redo the experiment again a decrease of the Height coefficient will increase the intercept</u>, and vice-versa, an increase of the Height coefficient is expected to drive down the intercept term (which is currently =-142.8215).
   plot(Height, Weight, main="Scatterplot Height vs. Weight", xlab="Height ", ylab="Weight ", pch=19)
@@ Line 237: / Line 237: @@
   coeffs
-    Estimate Std..Error        df   t.value     Pr...t..   df.Satt
+               Estimate  Std..Error   df     t.value     Pr...t..   df.Satt
   (Intercept) -143.26916 18.9698712  977.2059 -7.552458 9.792167e-14  977.2059
   Height         4.70294  0.2567663 1030.9565 18.316031   0.000000e+00 1030.9565
-               p.Satt
+                p.Satt
-  (Intercept) 9.792167e-14 Height      0.000000e+00
+  (Intercept)    9.792167e-14

@@ Line 199: / Line 199: @@
   Residual             295.652  17.195
   Number of obs: 1034, groups:  Team, 30; Position, 9
   Fixed effects:
   Estimate Std. Error t value
@@ Line 207: / Line 208: @@
     (Intr)
   Height -0.991
 The output starts with the user-specified model and data. Then, there’s the restricted maximum likelihood (ReML) criterion, summary of the model residuals (ε), followed by random and fixed effects.
@@ Line 291: / Line 242: @@
                 p.Satt
   (Intercept) 9.792167e-14 Height      0.000000e+00