0928hw2

1 load package

library(nlme)

2 input data

data(ergoStool, package="nlme")
str(ergoStool)

## 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(ergoStool)

## Grouped Data: effort ~ Type | Subject
##   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 descriptive statistc

# mean by Type
aggregate(ergoStool[,1], list(ergoStool$Type), mean)

##   Group.1         x
## 1      T1  8.555556
## 2      T2 12.444444
## 3      T3 10.777778
## 4      T4  9.222222

# mean by Subject
aggregate(ergoStool[,1], list(ergoStool$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

4 Linear mixed-effects model fit by REML

library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:nlme':
## 
##     lmList

summary(m0 <- lmer(effort ~ Type + (1 | Subject), data = ergoStool))

## Linear mixed model fit by REML ['lmerMod']
## Formula: effort ~ Type + (1 | Subject)
##    Data: ergoStool
## 
## REML criterion at convergence: 121.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.80200 -0.64317  0.05783  0.70100  1.63142 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1.775    1.332   
##  Residual             1.211    1.100   
## Number of obs: 36, groups:  Subject, 9
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   8.5556     0.5760  14.853
## TypeT2        3.8889     0.5187   7.498
## TypeT3        2.2222     0.5187   4.284
## TypeT4        0.6667     0.5187   1.285
## 
## Correlation of Fixed Effects:
##        (Intr) TypeT2 TypeT3
## TypeT2 -0.450              
## TypeT3 -0.450  0.500       
## TypeT4 -0.450  0.500  0.500

# fix effect 有 p 值的寫法
summary(m1 <- nlme::lme(effort ~ Type, random = ~ 1 | Subject, data=ergoStool, method="REML") )

## Linear mixed-effects model fit by REML
##  Data: ergoStool 
##        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

4.1 Approximate 95% confidence intervals

# 有fix and random coef est, 2.5%, 7.5%值
nlme::intervals(m1)

## 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

# 有fix and random 但無coef est值
confint(m0, method="boot")

## Computing bootstrap confidence intervals ...

## 
## 3 message(s): boundary (singular) fit: see ?isSingular

##                  2.5 %   97.5 %
## .sig01       0.4667008 2.135507
## .sigma       0.7535712 1.388692
## (Intercept)  7.4232996 9.714995
## TypeT2       2.8897877 4.935881
## TypeT3       1.1884255 3.313800
## TypeT4      -0.2838185 1.812693