1 Description

An experimenter wishes to compare 4 particular type of stools. Each of nine individuals randomly selected from a subject pool is asked to test the stools by recording the effort required to arise from each of 4 types of stools.

2 Input Data

data(ergoStool, package="nlme")
dta2 <- ergoStool
str(dta2)
## Classes 'nffGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   36 obs. of  3 variables:
##  $ effort : num  12 15 12 10 10 14 13 12 7 14 ...
##  $ Type   : Factor w/ 4 levels "T1","T2","T3",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ Subject: Ord.factor w/ 9 levels "8"<"5"<"4"<"9"<..: 8 8 8 8 9 9 9 9 6 6 ...
##  - attr(*, "formula")=Class 'formula'  language effort ~ Type | Subject
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Type of stool"
##   ..$ y: chr "Effort required to arise"
##  - attr(*, "units")=List of 1
##   ..$ y: chr "(Borg scale)"
##  - attr(*, "FUN")=function (x)  
##   ..- attr(*, "source")= chr "function (x) mean(x, na.rm = TRUE)"
##  - attr(*, "order.groups")= logi TRUE
head(dta2)
##   effort Type Subject
## 1     12   T1       1
## 2     15   T2       1
## 3     12   T3       1
## 4     10   T4       1
## 5     10   T1       2
## 6     14   T2       2

3 Data Analysis

3.1 Descriptive statistics

aggregate(dta2[,1], list(dta2$Type), mean)
##   Group.1         x
## 1      T1  8.555556
## 2      T2 12.444444
## 3      T3 10.777778
## 4      T4  9.222222
aggregate(dta2[,1], list(dta2$Subject), mean)
##   Group.1     x
## 1       8  8.25
## 2       5  8.50
## 3       4  9.25
## 4       9 10.00
## 5       6 10.25
## 6       3 10.75
## 7       7 10.75
## 8       1 12.25
## 9       2 12.25

3.2 Linear mixed model

# Random effects (random effect term: Subject)
summary(m0 <- nlme::lme(effort ~ Type, random = ~ 1 | Subject, data=dta2, method="REML") ) 
## Linear mixed-effects model fit by REML
##  Data: dta2 
##        AIC      BIC    logLik
##   133.1308 141.9252 -60.56539
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:    1.332465 1.100295
## 
## Fixed effects: effort ~ Type 
##                Value Std.Error DF   t-value p-value
## (Intercept) 8.555556 0.5760123 24 14.853079  0.0000
## TypeT2      3.888889 0.5186838 24  7.497610  0.0000
## TypeT3      2.222222 0.5186838 24  4.284348  0.0003
## TypeT4      0.666667 0.5186838 24  1.285304  0.2110
##  Correlation: 
##        (Intr) TypeT2 TypeT3
## TypeT2 -0.45               
## TypeT3 -0.45   0.50        
## TypeT4 -0.45   0.50   0.50 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.80200345 -0.64316591  0.05783115  0.70099706  1.63142054 
## 
## Number of Observations: 36
## Number of Groups: 9
faraway::sumary(m1 <-lme4::lmer(effort ~ Type + (1 | Subject), data=dta2) )
## Fixed Effects:
##             coef.est coef.se
## (Intercept) 8.56     0.58   
## TypeT2      3.89     0.52   
## TypeT3      2.22     0.52   
## TypeT4      0.67     0.52   
## 
## Random Effects:
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 1.33    
##  Residual             1.10    
## ---
## number of obs: 36, groups: Subject, 9
## AIC = 133.1, DIC = 123.2
## deviance = 122.1
faraway::sumary(m2 <-lme4::lmer(effort ~ 1 + (1 | Subject), data=dta2) )
## Fixed Effects:
## coef.est  coef.se 
##    10.25     0.48 
## 
## Random Effects:
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 1.03    
##  Residual             2.02    
## ---
## number of obs: 36, groups: Subject, 9
## AIC = 163.8, DIC = 158.5
## deviance = 158.2

3.3 Approximate 95% Confidence Intervals

nlme::intervals(m0)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                  lower      est.    upper
## (Intercept)  7.3667247 8.5555556 9.744386
## TypeT2       2.8183781 3.8888889 4.959400
## TypeT3       1.1517114 2.2222222 3.292733
## TypeT4      -0.4038442 0.6666667 1.737177
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Subject 
##                    lower     est.    upper
## sd((Intercept)) 0.749509 1.332465 2.368835
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.8292494 1.1002946 1.4599324
# fixed effects parameters only
confint(m1,parm="beta_", oldNames=F) 
## Computing profile confidence intervals ...
##                  2.5 %   97.5 %
## (Intercept)  7.4238425 9.687269
## TypeT2       2.8953043 4.882473
## TypeT3       1.2286377 3.215807
## TypeT4      -0.3269179 1.660251
# variance covariance
confint(m1,parm="theta_", oldNames=F)
## Computing profile confidence intervals ...
##                            2.5 %   97.5 %
## sd_(Intercept)|Subject 0.7342354 2.287261
## sigma                  0.8119798 1.390104
# plot
VCA::varPlot(effort ~ Type/Subject, 
             Data = dta2, 
             keep.order = T, 
             Points = list(pch=1, cex=1, col="gray"), 
             YLabel = list(text="effort",side=2,cex=1), 
             MeanLine = list(var=c("Type", "Subject"),
                             col=c("deepskyblue4","coral"), 
                             lwd=c(1,2)))
## Warning in min(x): min 中沒有無漏失的引數; 回傳 Inf
## Warning in max(x): max 中沒有無漏失的引數；回傳 -Inf