Drugs example

Eaach of 5 subjects is given 4 different drugs, with an appropriate “washout” time between drugs. The response is the reaction time to a stimulus.

##    Subject Drug Time
## 1        1    1   30
## 2        1    2   28
## 3        1    3   16
## 4        1    4   34
## 5        2    1   14
## 6        2    2   18
## 7        2    3   10
## 8        2    4   22
## 9        3    1   24
## 10       3    2   20
## 11       3    3   18
## 12       3    4   30
## 13       4    1   38
## 14       4    2   34
## 15       4    3   20
## 16       4    4   44
## 17       5    1   26
## 18       5    2   28
## 19       5    3   14
## 20       5    4   30

Fixed effects ANOVA

    par(mfrow = c(1., 2.),pty="s")
  boxcox(lm(Time~factor(Drug)+factor(Subject)))

  react1.ran <- aov(log(Time) ~ factor(Drug) + factor(Subject))
    print(summary(react1.ran))

##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Drug)     3  1.395  0.4650   35.08 3.22e-06 ***
## factor(Subject)  4  1.206  0.3015   22.75 1.58e-05 ***
## Residuals       12  0.159  0.0133                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  par(mfrow = c(1., 2.),pty="s")
    plot(Subject, log(Time))
    points(Subject, react1.ran$fitted, pch = "+",col="blue")
    plot(Drug, log(Time))
    points(Drug, react1.ran$fitted, pch = "+",col="blue")

Residuals: \(e_{ij}=y_{ij}-\hat{y}_{ij}=y_{ij}-\hat{\mu}-\hat{\alpha}_i-\hat{\beta}_j=y_{ij}-(\bar{y}_{i .}-\bar{y})-(\bar{y}_{. j}-\bar{y})=y_{ij}-\bar{y}_{i .}-\bar{y}_{. j}+\bar{y}\)

Mixed effects ANOVA

#   react2.ran <- aov(log(Time) ~ factor(Drug) + Error(factor(Subject)))
    react2.ran <- lme(log(Time) ~ factor(Drug), random=~1|Subject)
    print(summary(react2.ran,correlation=F))

## Linear mixed-effects model fit by REML
##  Data: NULL 
##        AIC      BIC   logLik
##   7.165977 11.80151 2.417011
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:   0.2684461 0.1151237
## 
## Fixed effects: log(Time) ~ factor(Drug) 
##                   Value  Std.Error DF   t-value p-value
## (Intercept)    3.222798 0.13062678 12 24.671803  0.0000
## factor(Drug)2 -0.007424 0.07281062 12 -0.101957  0.9205
## factor(Drug)3 -0.502731 0.07281062 12 -6.904642  0.0000
## factor(Drug)4  0.217999 0.07281062 12  2.994058  0.0112
##  Correlation: 
##               (Intr) fc(D)2 fc(D)3
## factor(Drug)2 -0.279              
## factor(Drug)3 -0.279  0.500       
## factor(Drug)4 -0.279  0.500  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.63034180 -0.40049189 -0.07857529  0.49359790  1.75685990 
## 
## Number of Observations: 20
## Number of Groups: 5

    print(anova(react2.ran))

##              numDF denDF  F-value p-value
## (Intercept)      1    12 658.0937  <.0001
## factor(Drug)     3    12  35.0825  <.0001

    par(mfrow = c(1., 2.),pty="s")
    plot(Drug,log(Time))
    points(Drug,react2.ran$fitted[,1],pch="+",col="blue")
      plot(fitted(react2.ran), resid(react2.ran))

Residuals: \(e_{ij}=y_{ij}-\hat{y}_{ij}=y_{ij}-\hat{\mu}-\hat{\alpha}_i-\hat{\beta}_j=y_{ij}-\bar{y}_{. j}-\frac{\sigma^2_\alpha}{\sigma^2_\alpha+\sigma^2/n_i}(\bar{y}_{i .}-\bar{y})\)

Random effects: \(e_i^\alpha=\hat{\alpha}_i\)

par(mfrow = c(1., 2.),pty="s")
      qqnorm(react2.ran$res,main = "Normal Q-Q plot for residuals")
      qqnorm(ranef(react2.ran)[,1],main="Normal Q-Q plot for random effects")