#----------------------------------------------------------------------------------------
#Statistical Computing Seminars Repeated Measures Analysis using SAS
#http://www.ats.ucla.edu/stat/sas/seminars/sas_repeatedmeasures/default.htm
#----------------------------------------------------------------------------------------

################################################
#		Demo Analysis #1
#被験者内1因子，被験者間1因子（3 * 2）の分散分析
################################################
##(a)データの作成
dat1 <- matrix(c(
  1,        10,       10,       10,
  1,        10,       10,       10,
  1,        10,       10,       10,
  1,        10,       10,       10,
  2,        15,       15,       15,
  2,        15,       15,       15,
  2,        16,       15,       15,
  2,        15,       15,       15
  ),byrow=T,ncol=4)

dat1 <- data.frame(dat1)
colnames(dat1) <- c("group", "time1", "time2", "time3")
rownames(dat1) <- paste("s",seq(1:nrow(dat1)),sep="")
dat1$group <- factor(dat1$group)


##(b) 縦断データ用に整形
dat1s <- stack(dat1)
dat1s$group <- rep(dat1$group,3)
dat1s$subject <- factor(rep(rownames(dat1),3))
names(dat1s) <- c("values","time","group","subject")


##(c) 基礎データ
attach(dat1s)
(fig1 <- tapply(values,list(group,time),mean))	#平均値
tapply(values,list(group,time),length)		#標本サイズ
matplot(t(fig1),pch=c(2,4),cex=1.4,xlab="Time",ylab="Means")
matlines(t(fig1),lwd=2,lty=1)

graphics.off()					#図を閉じる
detach(dat1s)


##(d) 分散分析
aov.dat1 <- aov(values~ group * time +  Error(subject/time) ,data=dat1s)
summary(aov.dat1)



################################################
#		Demo Analysis #2
#被験者内1因子，被験者間1因子（3 * 2）の分散分析
################################################

##(a)データの作成
dat2 <- matrix(c(
  1,        14,       19,       29,
  1,        15,       25,       26,
  1,        16,       16,       31,
  1,        12,       24,       32,
  2,        10,       21,       24,
  2,        17,       26,       35,
  2,        19,       22,       32,
  2,        15,       23,       34
),byrow=T,ncol=4)
dat2 <- data.frame(dat2)
colnames(dat2) <- c("group", "time1", "time2", "time3")
rownames(dat2) <- paste("s",seq(1:nrow(dat2)),sep="")
dat2$group <- factor(dat2$group)


##(b) 縦断データ用に整形
dat2s <- stack(dat2)
dat2s$group <- rep(dat2$group,3)
dat2s$subject <- factor(rep(rownames(dat2),3))
names(dat2s) <- c("values","time","group","subject")


##(c) 基礎データ
attach(dat2s)
(fig1 <- tapply(values,list(group,time),mean))	#平均値
tapply(values,list(group,time),length)		#標本サイズ
matplot(t(fig1),pch=c(2,4),cex=1.4,xlab="Time",ylab="Means")
matlines(t(fig1),lwd=2,lty=1)
graphics.off()
detach(dat2s)


##(d) 分散分析
aov.dat2 <- aov(values~ group * time +  Error(subject/time) ,data=dat2s)
summary(aov.dat2)



################################################
#		Demo Analysis #3
#被験者内1因子，被験者間1因子（3 * 2）の分散分析
################################################

##(a)データの作成
dat3 <- matrix(c(
  1,        35,       25,       16,       
  1,        32,       23,       12,       
  1,        36,       22,       14,       
  1,        34,       21,       13,       
  2,        57,       43,       22,       
  2,        54,       46,       26,       
  2,        55,       46,       23,       
  2,        60,       47,       25       
),byrow=T,ncol=4)
dat3 <- data.frame(dat3)
colnames(dat3) <- c("group", "time1", "time2", "time3")
rownames(dat3) <- paste("s",seq(1:nrow(dat3)),sep="")
dat3$group <- factor(dat3$group)


##(b) 縦断データ用に整形
dat3s <- stack(dat3)
dat3s$group <- rep(dat3$group,3)
dat3s$subject <- factor(rep(rownames(dat3),3))
names(dat3s) <- c("values","time","group","subject")


##(c) 基礎データ
attach(dat3s)
(fig1 <- tapply(values,list(group,time),mean))	#平均値
tapply(values,list(group,time),length)		#標本サイズ
matplot(t(fig1),pch=c(2,4),cex=1.4,xlab="Time",ylab="Means")
matlines(t(fig1),lwd=2,lty=1)
graphics.off()
detach(dat3s)


##(d) 分散分析
aov.dat3 <- aov(values~ group * time +  Error(subject/time) ,data=dat3s)
summary(aov.dat3)



################################################
#		Demo Analysis #4
#被験者内1因子，被験者間1因子（3 * 2）の分散分析
################################################

##(a)データの作成
dat4 <- matrix(c(
  1,        35,       25,       12,
  1,        34,       22,       13,
  1,        36,       21,       18,
  1,        35,       23,       15,
  2,        31,       43,       57,
  2,        35,       46,       58,
  2,        37,       48,       51,
  2,        32,       45,       53
),byrow=T,ncol=4)
dat4 <- data.frame(dat4)
colnames(dat4) <- c("group", "time1", "time2", "time3")
rownames(dat4) <- paste("s",seq(1:nrow(dat4)),sep="")
dat4$group <- factor(dat4$group)


##(b) 縦断データ用に整形
dat4s <- stack(dat4)
dat4s$group <- rep(dat4$group,3)
dat4s$subject <- factor(rep(rownames(dat4),3))
names(dat4s) <- c("values","time","group","subject")


##(c) 基礎データ
attach(dat4s)
(fig1 <- tapply(values,list(group,time),mean))	#平均値
tapply(values,list(group,time),length)		#標本サイズ
matplot(t(fig1),pch=c(2,4),cex=1.4,xlab="Time",ylab="Means")
matlines(t(fig1),lwd=2,lty=1)
graphics.off()
detach(dat4s)


##(d) 分散分析
aov.dat4 <- aov(values~ group * time +  Error(subject/time) ,data=dat4s)
summary(aov.dat4)



################################################
#	Exercise data examples
#被験者内1因子，被験者間2因子（3 * 2 * 3）の分散分析
#
#	 uncomplete job: manova, xlab
################################################

##(a)データの作成
dat5 <- matrix(c(
1,     1,         1,       85,       85,       88,
2,     1,         1,       90,       92,       93,
3,     1,         1,       97,       97,       94,
4,     1,         1,       80,       82,       83,
5,     1,         1,       91,       92,       91,
6,     1,         2,       83,       83,       84,
7,     1,         2,       87,       88,       90,
8,     1,         2,       92,       94,       95,
9,     1,         2,       97,       99,       96,
10,    1,         2,      100,       97,      100,
11,    2,         1,       86,       86,       84,
12,    2,         1,       93,      103,      104,
13,    2,         1,       90,       92,       93,
14,    2,         1,       95,       96,      100,
15,    2,         1,       89,       96,       95,
16,    2,         2,       84,       86,       89,
17,    2,         2,      103,      109,       90,
18,    2,         2,       92,       96,      101,
19,    2,         2,       97,       98,      100,
20,    2,         2,      102,      104,      103,
21,    3,         1,       93,       98,      110,
22,    3,         1,       98,      104,      112,
23,    3,         1,       98,      105,       99,
24,    3,         1,       87,      132,      120,
25,    3,         1,       94,      110,      116,
26,    3,         2,       95,      126,      143,
27,    3,         2,      100,      126,      140,
28,    3,         2,      103,      124,      140,
29,    3,         2,       94,      135,      130,
30,    3,         2,       99,      111,      150
),byrow=T,ncol=6)

dat5 <- data.frame(dat5)
colnames(dat5) <- c("id", "exertype","diet", "time1", "time2", "time3")
rownames(dat5) <- paste("s",seq(1:nrow(dat5)),sep="")



##(b) 縦断データ用に整形
dat5s <- stack(dat5[,c(-1,-2,-3)])
dat5s$exertype <- factor(rep(dat5$exertype,3))
dat5s$diet <- factor(rep(dat5$diet,3))
dat5s$subject <- factor(rep(rownames(dat5),3))
names(dat5s) <- c("values","time","exertype","diet","subject")


##(c) 基礎データ
attach(dat5s)
(fig1 <- tapply(values,list(diet,time),mean))		#平均値
(fig2 <- tapply(values,list(exertype,time),mean))	#平均値
tapply(values,list(diet,time),length)			#標本サイズ
tapply(values,list(exertype,time),length)		#標本サイズ

matplot(t(fig1),pch=c(2,4),cex=1.4,xlab="Time",ylab="Means")
matlines(t(fig1),lwd=2,lty=1)
graphics.off()

matplot(t(fig2),pch=c(2,4),cex=1.4,xlab="Time",ylab="Means")
matlines(t(fig2),lwd=2,lty=1)
graphics.off()

detach(dat5s)


##(d) 分散分析 (model1)
aov.dat5 <- aov(values~ diet * time +  Error(subject/time) ,data=dat5s)
summary(aov.dat5)


##(e) 分散分析 (model2)
aov.dat5_2 <- aov(values~ exertype * time +  Error(subject/time) ,data=dat5s)
summary(aov.dat5_2)