# 
# c-stat example
#

# read in data 
dta <- read.table("cstat.txt", header=T)

str(dta)
## 'data.frame':    42 obs. of  1 variable:
##  $ nc: int  28 46 39 45 24 20 35 37 36 40 ...
head(dta)
##   nc
## 1 28
## 2 46
## 3 39
## 4 45
## 5 24
## 6 20
#
# plot figure 1
#
plot(1:42, dta[,1], xlab="Observations", ylab="Number of Children")
lines(1:42, dta[,1])
abline(v=10, lty=2)
abline(v=32, lty=2)
segments(1, mean(dta[1:10,1]),10, mean(dta[1:10,1]),col="red")
segments(11, mean(dta[11:32,1]),32, mean(dta[11:32,1]),col="red")
segments(33, mean(dta[33:42,1]),42, mean(dta[33:42,1]),col="red")

#
# calculate c-stat for first baseline phase
# 帶入function指令

dat <- function(x){
 cden <- 1-(sum(diff(x)^2)/(2*(10-1)*var(x)))
 sc <- sqrt((10-2)/((10-1)*(10+1)))
pval <- 1-pnorm(cden/sc)
return(pval)}

dat(dta[1:10,1])
## [1] 0.2866238
#
# calculate c-stat for first baseline plus group tokens
# 帶入function指令

data <- function(y){
  n <- 32
  cden <- 1-(sum(diff(y)^2)/(2*(n-1)*var(y)))
  sc <- sqrt((n-2)/((n-1)*(n+1)))
  pval <- 1-pnorm(cden/sc)
  list <- list(z=cden/sc,pvalue=pval)
return(list)}

data(dta[1:32,1])
## $z
## [1] 3.879054
## 
## $pvalue
## [1] 5.243167e-05