STAT 700, Fall 2014 Homework 6 Problems due Wed. Nov. 12

• Francis Alvarez
• Gonzalo Urritia

require(nlme)

## Loading required package: nlme

Ortho=Orthodont
OrthoM <- Orthodont[Orthodont$Sex == "Male",] #Only include male patients
par(pty="s")

1. Fit two mixed models, fit1 and fit2, one with a random effect for the intercept and one with random effects for both the intercept and slope. Provide the summary for each model.

fit1 <- lme(distance~age, data=OrthoM, random=~1 | Subject) #Strictly random effect for intercept
summary(fit1)

## Linear mixed-effects model fit by REML
##  Data: OrthoM 
##     AIC BIC logLik
##   281.4 290 -136.7
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:       1.625    1.678
## 
## Fixed effects: distance ~ age 
##              Value Std.Error DF t-value p-value
## (Intercept) 16.341    1.1287 47  14.477       0
## age          0.784    0.0938 47   8.361       0
##  Correlation: 
##     (Intr)
## age -0.914
## 
## Standardized Within-Group Residuals:
##       Min        Q1       Med        Q3       Max 
## -3.008055 -0.640689  0.007833  0.534481  3.052947 
## 
## Number of Observations: 64
## Number of Groups: 16

fit2 <- lme(distance~age, data=OrthoM, random=~age | Subject) #If wanted just random effect for slope would do random=~age-1
summary(fit2)

## Linear mixed-effects model fit by REML
##  Data: OrthoM 
##     AIC   BIC logLik
##   285.1 297.9 -136.6
## 
## Random effects:
##  Formula: ~age | Subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev Corr  
## (Intercept) 2.6666 (Intr)
## age         0.1887 -0.788
## Residual    1.6091       
## 
## Fixed effects: distance ~ age 
##              Value Std.Error DF t-value p-value
## (Intercept) 16.341    1.2099 47  13.506       0
## age          0.784    0.1016 47   7.722       0
##  Correlation: 
##     (Intr)
## age -0.926
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -2.65543 -0.59382  0.01946  0.57552  3.15772 
## 
## Number of Observations: 64
## Number of Groups: 16

2. Use the two models fit1 and fit2 test \(H_0: \sigma^{2}_{slope} \ne {0}\) against \(H_1: \sigma^{2}_{slope} = {0}\). State your conclusion and state which is the “best” model.

anova(fit2, fit1)

##      Model df   AIC   BIC logLik   Test L.Ratio p-value
## fit2     1  6 285.1 297.9 -136.6                       
## fit1     2  4 281.4 290.0 -136.7 1 vs 2  0.3106  0.8561

After looking at the anova test we see that we fail to reject the null hypothesis (due to such a high p-value=0.8561) and accept the fit1 as the better model. We can also look at the AIC and BIC, which shows that fit1 is lower suggesting that it is the better model.

3. Using your answer to part (b). Give the estimates for all parameters in the “best” model.

The parameters in the “best model” are:

Fixed Intercept: \(\hat{\beta_0}=16.341\) Fixed Age(slope): \(\hat{\beta_0}=0.784\) Random Effect Intercept: \(\hat{\sigma}^{2}_{intercept}= (1.625)^2 = 2.6406\) Residual: \(\hat{\sigma}^{2}_{intercept}= (0.7800)^2 = 0.6084\)

4. How well does the “best” model fit the data? Include and examine the diagnostics plot of the residuals. You can also look at a Q-Q plot of the residuals.

plot(fit1)

The residual plot is scattered and does not have patterns so the model is a good fit.

5. Use the R augPred function to plot the fitted values for the “best” model.

plot(augPred(fit1))

The model fits the data extremely well. The only exception is Male09 but his data differs greatly from the rest of the group so in that sense it is understable. It may be in the best interest to remove him from the data set, otherwise the model is a good fit.