ANOVAQ1-3
Li
2017年10月22日
Q1
#random-effects analysis
require(pacman)## Loading required package: pacmanpacman::p_load(HLMdiag,lme4,lattice)
#lme4 for lmer
dta <- radon
head(dta)## log.radon basement uranium county county.name
## 1 0.7884574 1 -0.6890476 1 AITKIN
## 2 0.7884574 0 -0.6890476 1 AITKIN
## 3 1.0647107 0 -0.6890476 1 AITKIN
## 4 0.0000000 0 -0.6890476 1 AITKIN
## 5 1.1314021 0 -0.8473129 2 ANOKA
## 6 0.9162907 0 -0.8473129 2 ANOKAstr(dta)## 'data.frame': 919 obs. of 5 variables:
## $ log.radon : num 0.788 0.788 1.065 0 1.131 ...
## $ basement : int 1 0 0 0 0 0 0 0 0 0 ...
## $ uranium : num -0.689 -0.689 -0.689 -0.689 -0.847 ...
## $ county : int 1 1 1 1 2 2 2 2 2 2 ...
## $ county.name: Factor w/ 85 levels "AITKIN","ANOKA",..: 1 1 1 1 2 2 2 2 2 2 ...#random-effects analysis of variance on radon by county
summary(m1 <- lmer(log.radon ~ (1|county.name) , data = dta))## Linear mixed model fit by REML ['lmerMod']
## Formula: log.radon ~ (1 | county.name)
## Data: dta
##
## REML criterion at convergence: 2259.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4661 -0.5734 0.0441 0.6432 3.3516
##
## Random effects:
## Groups Name Variance Std.Dev.
## county.name (Intercept) 0.09581 0.3095
## Residual 0.63662 0.7979
## Number of obs: 919, groups: county.name, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.31258 0.04891 26.84str(fitted(m1))## Named num [1:919] 1.067 1.067 1.067 1.067 0.888 ...
## - attr(*, "names")= chr [1:919] "1" "2" "3" "4" ...histogram(~ fitted(m1), type = "density", nint=23,
xlim = c(0.5, 2), ylim = c(0, 8),
xlab ="Fitted value", main = "radom-effects")
Q2
#mixed effects model
require(nlme)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmListdta<-ergoStool
str(dta)## 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#mixed effects model
summary(me <- lme(effort ~ Type ,random = ~ 1 | Subject, data = dta))## Linear mixed-effects model fit by REML
## Data: dta
## 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: 9intervals(me)## 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.4599324Q3
#One Way Repeated measures
setwd("C:\\Users\\USER\\Desktop\\Multilevel Statistical\\w5 exercise")
dta <- read.table(file='thetaEEG.txt',header=TRUE)
str(dta)## 'data.frame': 19 obs. of 10 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Ch3 : num -3.54 5.72 0.52 0 2.07 1.67 9.13 -0.43 -0.56 1.28 ...
## $ Ch4 : num -3.11 5.07 -0.18 0.74 0.76 ...
## $ Ch5 : num -0.24 6.87 0.9 1.1 3.51 2.77 3.44 -0.31 -1.22 1.89 ...
## $ Ch6 : num 0.42 5.96 0.6 0.13 0.6 4.55 4.8 -0.61 0.67 1.77 ...
## $ Ch7 : num -0.49 8.2 1.27 0.19 3.71 1.8 0.48 -1.04 -0.97 1.83 ...
## $ Ch8 : num 2.13 4.87 1.28 0.07 1.86 4.79 1.63 -0.13 -0.98 0.91 ...
## $ Ch17: num -4.15 5.48 -0.95 0.8 1.49 4.51 9.94 -0.61 0 1.4 ...
## $ Ch18: num 2.87 5.57 1.74 0.25 3.11 3.24 1.34 -0.61 -1.22 1.1 ...
## $ Ch19: num 1.34 6.33 0.79 -0.66 1.8 3.99 1.53 -0.43 -0.91 -0.12 ...#set the digits
options(digits=3)
#channel mean
colMeans(dta[,2:10])## Ch3 Ch4 Ch5 Ch6 Ch7 Ch8 Ch17 Ch18 Ch19
## 0.901 1.589 1.037 1.146 0.851 0.853 1.403 0.854 0.995#cov
cov(dta[,2:10])## Ch3 Ch4 Ch5 Ch6 Ch7 Ch8 Ch17 Ch18 Ch19
## Ch3 8.65 9.53 5.63 5.21 3.85 2.91 9.02 3.55 3.62
## Ch4 9.53 12.35 5.21 6.17 2.67 3.22 10.72 2.42 4.04
## Ch5 5.63 5.21 5.56 4.21 4.54 3.50 5.73 4.22 4.24
## Ch6 5.21 6.17 4.21 5.07 3.08 3.71 6.07 3.66 4.14
## Ch7 3.85 2.67 4.54 3.08 5.18 2.90 3.45 3.81 3.24
## Ch8 2.91 3.22 3.50 3.71 2.90 4.28 3.59 3.19 3.82
## Ch17 9.02 10.72 5.73 6.07 3.45 3.59 10.74 3.38 4.46
## Ch18 3.55 2.42 4.22 3.66 3.81 3.19 3.38 6.77 2.64
## Ch19 3.62 4.04 4.24 4.14 3.24 3.82 4.46 2.64 5.57#reshape the data
require(reshape2)## Loading required package: reshape2#measure.vars=要併的變項範圍 var.ids=併的變項名稱
#variable.name= 併成的變項 value.name=併成的變項值名稱
d<-melt(dta, measure.vars= c(2:10),var.ids=c("ch3", "ch4", "ch5","ch6","ch7" , "ch8","ch17","ch18","ch19"),
variable.name="ch", value.name="sc")
#require(tidyr) # d<- gather(dta,key = Channel, value =sc ,2:10)
library(car)
#One Way Repeated measures Anova
summary(aov(sc~ ch + Error(ID / ch), d))##
## Error: ID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 112 112
##
## Error: ID:ch
## Df Sum Sq Mean Sq
## ch 8 23.1 2.89
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## ch 8 10 1.23 0.18 0.99
## Residuals 153 1020 6.67