library(ggplot2)
library(lme4)
library(dplyr)

#Dataset Creation:

#Pre/Post Dataset:
  #New alarm
before_1 <- round(rnorm(n=20, mean=20, sd=3),1)
diff_1   <- round(rnorm(n=20, mean=-2, sd=3),1) 
after_1  <- before_1 + diff_1

print(before_1); print(diff_1); print(after_1)

  #Control
before_0 <- round(rnorm(n=20, mean=20, sd=3),1)
diff_0   <- round(rnorm(n=20, mean=0, sd=3),1) 
after_0  <- before_0 + diff_0

print(before_0); print(diff_0); print(after_0)

alarm_df1 <-data.frame('ID'    =c(rep(seq(1,20),2), rep(seq(21,40),2)),
                      'alarm' =c(rep('Control',40),rep('New_Alarm',40)),
                      'time'  =c(rep('Before',20), rep('Post',20), rep('Before',20), rep('Post',20)),
                      'snooze'=c(before_0, after_0, before_1, after_1))


alarm_df2 <-data.frame('ID'    =c(seq(1,40)),
                       'alarm' =c(rep('Control',20),rep('New_Alarm',20)),
                       'snooze_pre'=c(before_0, before_1),
                       'snooze_post'=c(after_0, after_1),
                       'snooze_diff'=c(after_0-before_0, after_1-before_1))


#Repeated Measures ANOVA:
sleep_0 <-round(rnorm(n=60, mean=7, sd=1),1) #control
sleep_1 <-c(round(rnorm(n=20, mean=7, sd=1),1), round(rnorm(n=20, mean=7.5, sd=1),1), 
            round(rnorm(n=20, mean=8, sd=1),1)) #mattress
ID <-c(rep(seq(1,20),3), rep(seq(21,40),3))
month <-c(rep('1',20), rep('2',20), rep('3',20), rep('1',20), rep('2',20), rep('3',20))
treat <-c(rep('control',60), rep('mattress',60))

bed_df1 <-data.frame('ID'=ID,
                     'sleep'=c(sleep_0, sleep_1),
                     'month'=month,
                     'treat'=treat)

#Mixed Models Regression:
ID    <-c(rep(seq(1,10),5))
month <-c(rep(1,10), rep(2,10), rep(3,10), rep(4,10), rep(5,10))
brand <-c(rep(c('A','A','A','A','A','B','B','B','B','B'),5))
sleep <-c(round(rnorm(n=10, mean=7.0, sd=1),1), round(rnorm(n=10, mean=7.3, sd=1),1), 
          round(rnorm(n=10, mean=7.5, sd=1),1), round(rnorm(n=10, mean=7.8, sd=1),1),
          round(rnorm(n=10, mean=8.0, sd=1),1))

bed_df2 <-data.frame('ID'=ID,
                    'sleep'=sleep,
                    'month'=month,
                    'brand'=brand)

#Data already made
#alarm_df1 <- read.csv("C:\\...\\alarm_df1.csv")
#alarm_df2 <- read.csv( "C:\\...\\2022 Data\\alarm_df2.csv")
#bed_df1   <- read.csv("C:\\...\\bed_df1.csv")
#bed_df2   <- read.csv("C:\\...\\bed_df2.csv")


#########################################################################################################

#Exercise 1: Pre/post test

  #1.1: no control
head(alarm_df1)
alarm_df_sub <-alarm_df1[alarm_df1$alarm=="New_Alarm",]

boxplot(data=alarm_df_sub, snooze~time)
t.test(data=alarm_df_sub, snooze~time, paired=T)
ggplot(data=alarm_df_sub, aes(x=time, y=snooze, group=ID)) + geom_line()

  #1.2: control
head(alarm_df2)
boxplot(data=alarm_df2, snooze_diff~alarm)
t.test(alarm_df2$snooze_diff[alarm_df2$alarm=='Control'],
       alarm_df2$snooze_diff[alarm_df2$alarm=='New_Alarm'])
ggplot(data=alarm_df1, aes(x=time, y=snooze, group=ID, color=alarm)) + geom_line()


#Exercise 2: Repeated-Measures ANOVA

  #2.1: Month alone
head(bed_df1)
bed_df1$month=as.factor(bed_df1$month)
boxplot(data=bed_df1, sleep~month)

aov1 <-lm(data=bed_df1, sleep~month)
summary(aov1)

aov2 <-lmer(data=bed_df1, sleep~month + (1|ID))
summary(aov2)

bed_sum<-bed_df1 %>% 
  group_by(month) %>%       
  summarise(mean=mean(sleep),
            sd = sd(sleep),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
ggplot(data=bed_sum, aes(x=month, y=mean, fill=month)) +
  geom_bar(stat="identity") +
  scale_y_continuous(limits =c(0,8.5)) +
  geom_errorbar(aes(ymin=ll, ymax=ul), width=0.1)



  #2.2 Month and Treatment
boxplot(data=bed_df1, sleep~treat*month, col=c("red","blue"))

aov3 <-lm(data=bed_df1, sleep~month*treat)
summary(aov3)

aov4 <-lmer(data=bed_df1, sleep~month*treat + (1|ID))
summary(aov4)

bed_sum2<-bed_df1 %>% 
  group_by(month, treat) %>%      
  summarise(mean=mean(sleep),
            sd = sd(sleep),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
ggplot(data=bed_sum2, aes(x=month, y=mean, fill=treat)) +
  geom_bar(aes(fill=treat), stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ll, ymax=ul, group=treat), width=0.1, position=position_dodge(0.9)) +
  geom_text(aes(label="*"), x=c(3.225), y=c(8.80), size=7) 


#Exercise 3: Mixed Models

  #3.1 Month alone
head(bed_df2)
plot(bed_df2$month, bed_df2$sleep)

mm1 <-lm(data=bed_df2, sleep~month)
summary(mm1)

mm2 <-lmer(data=bed_df2, sleep~month + (1|ID))
summary(mm2)

ggplot(data=bed_df2, aes(x=month, y=sleep)) +
  geom_point() + geom_smooth(method=lm, color="red")

  #3.2 Month and Brand
mm3 <-lm(data=bed_df2, sleep~month*brand)
summary(mm3)

mm4 <-lmer(data=bed_df2, sleep~month*brand + (1|ID))
summary(mm4)


#Output to SAS:
write.csv(alarm_df1, "C:\\...\\alarm_df1.csv")
write.csv(alarm_df2, "C:\\...\\alarm_df2.csv")
write.csv(bed_df1, "C:\\...\\2022 Data\\bed_df1.csv")
write.csv(bed_df2, "C:\\...\\bed_df2.csv")



###################

#Y-variable data replicated
#alarm_df1$snooze
sn1<-c(24.0, 21.3, 23.5, 21.7, 19.2, 24.3, 19.3, 16.5, 18.9, 16.7, 21.7, 23.6, 25.8, 17.7, 18.9, 24.2, 20.2, 22.2, 23.1, 22.0, 
       18.2, 17.7, 24.1, 25.9, 16.6, 23.5, 22.0, 15.5, 21.2, 17.3, 22.1, 22.5, 19.8, 18.2, 18.7, 22.8, 22.0, 20.9, 19.9, 22.0, 
       23.0, 22.3, 18.1, 18.6, 18.4, 18.6, 25.9, 15.8, 22.8, 16.6, 17.6, 23.3, 22.6, 17.7, 18.9, 21.1, 19.7, 21.1, 20.9, 17.7, 
       23.1, 22.4, 14.9, 17.0, 20.7, 19.0, 22.7, 16.0, 23.8, 13.4, 18.0, 24.4, 19.4, 19.8, 15.3, 19.3, 18.3, 18.2, 16.1, 14.4)

#alarm_df2$snooze_pre
spr<-c(24.0, 21.3, 23.5, 21.7, 19.2, 24.3, 19.3, 16.5, 18.9, 16.7, 21.7, 23.6, 25.8, 17.7, 18.9, 24.2, 20.2, 22.2, 23.1, 22.0, 
       23.0, 22.3, 18.1, 18.6, 18.4, 18.6, 25.9, 15.8, 22.8, 16.6, 17.6, 23.3, 22.6, 17.7, 18.9, 21.1, 19.7, 21.1, 20.9, 17.7)

#alarm_df2$snooze_post
spo<-c(18.2, 17.7, 24.1, 25.9, 16.6, 23.5, 22.0, 15.5, 21.2, 17.3, 22.1, 22.5, 19.8, 18.2, 18.7, 22.8, 22.0, 20.9, 19.9, 22.0, 
       23.1, 22.4, 14.9, 17.0, 20.7, 19.0, 22.7, 16.0, 23.8, 13.4, 18.0, 24.4, 19.4, 19.8, 15.3, 19.3, 18.3, 18.2, 16.1, 14.4)

#alarm_df2$snooze_diff
spd<-c(-5.8, -3.6, 0.6, 4.2, -2.6, -0.8, 2.7, -1.0, 2.3, 0.6, 0.4, -1.1, -6.0, 0.5, -0.2, -1.4, 1.8, -1.3, -3.2, 0.0, 0.1, 0.1,
       -3.2, -1.6, 2.3, 0.4, -3.2, 0.2, 1.0, -3.2, 0.4, 1.1, -3.2, 2.1, -3.6, -1.8, -1.4, -2.9, -4.8, -3.3)

#bed_df1$sleep
sl1<-c(7.4, 8.6, 6.1, 7.9, 8.0, 6.5, 7.1, 6.4, 6.4, 6.7, 7.1, 7.6, 8.0, 5.3, 8.4, 7.5, 8.0, 5.5, 7.2, 6.6, 7.9, 8.0, 5.2, 7.6, 6.2, 
       8.2, 8.2, 8.6, 6.1, 7.2, 8.3, 5.9, 8.2, 7.0, 7.1, 9.3, 6.6, 8.4, 6.6, 9.0, 7.1, 6.5, 6.1, 7.1, 6.1, 7.8, 6.2, 7.2, 5.7, 6.4, 
       7.0, 8.9, 9.7, 7.4, 6.5, 6.5, 5.8, 7.9, 6.9, 5.5, 6.6, 7.0, 7.4, 6.3, 6.7, 8.9, 6.0, 7.7, 5.1, 7.6, 5.8, 5.3, 6.1, 6.1, 7.2, 
       6.1, 6.3, 9.0, 6.4, 5.9, 9.4, 7.1, 8.6, 6.5, 7.7, 7.2, 7.8, 7.9, 5.3, 6.9, 6.5, 5.4, 6.1, 7.2, 8.3, 6.5, 7.9, 8.1, 6.9, 8.1, 
       7.6, 6.5, 7.8, 6.2, 7.1, 8.8, 7.3, 8.9, 8.9, 7.2, 7.5, 8.3, 7.5, 8.8, 7.6, 8.3, 7.0, 9.1, 7.3, 8.4)

#bed_df2$sleep
sl2<-c(6.0, 7.9, 7.7, 6.7, 6.2, 7.8, 6.5, 8.2, 7.6, 6.2, 7.9, 6.5, 8.2, 7.8, 6.4, 7.2, 6.0, 6.8, 8.5, 6.9, 9.2, 7.8, 7.3, 6.2, 6.6, 
       7.9, 6.4, 7.8, 7.4, 8.0, 9.5, 7.7, 8.4, 7.6, 7.9, 7.6, 7.9, 7.1, 7.4, 9.3, 8.7, 8.6, 8.7, 7.0, 8.3, 7.6, 8.1, 10.7, 5.8, 9.1)