setwd("~/Desktop")
library(ggplot2)
library(broom)
cls=c(lat="numeric",lon="numeric")
wt=read.table(text="255.9194
216.9117
277.9368
313.241
211.0513
254.0302
249.1868
229.3152
282.6375
200.3201
246.521
250.9596
282.7773
285.2515
275.1814
114.234
296.2902
211.8189
248.252
170.2063", stringsAsFactors=F,na.strings="unknown", header=F, sep="\t")
wt=as.numeric(unlist(wt))
sd_wt=sd(wt) #sd=46.99014
dmean1=0.4
dmean2=0.2
dmean3=0.1
thesis=0.75
mean_wt=mean(wt) #mean=243.6021
df=data.frame(group=rep("wild-type",20),action_potential_duration=wt)
hist(wt,xlab="action potential duration",main=NULL)
dev.copy(png,"Q4_wt_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
ggplot(df,aes(x=group,y=action_potential_duration,color=group))+labs(title=NULL,x=NULL,y="action potential duration (ms)",color=NULL) +geom_point(size=3)+ geom_point(stat="summary",color="black",size=4, shape=18)+geom_errorbar(aes(ymin=mean_wt-sd_wt,ymax=mean_wt+sd_wt),length=5,width=0.2)+theme(axis.text.x = element_blank(),text=element_text(size=24,face="bold"))
## Warning: Ignoring unknown parameters: length
## No summary function supplied, defaulting to `mean_se()`
dev.copy(png,"Q4_wt_data.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
shapiro.test(wt) #p-value=0.13
##
## Shapiro-Wilk normality test
##
## data: wt
## W = 0.92611, p-value = 0.13
#40% change in mean action potential duration
d1=(dmean1*mean_wt)/sd_wt
power40=tidy(power.t.test(n=NULL,sig.level = .05,power=0.8,sd=1,delta=d1))
#20% change in mean action potential duration
d2=(dmean2*mean_wt)/sd_wt
power20=tidy(power.t.test(n=NULL,sig.level = .05,power=0.8,sd=1,delta=d2))
#10% change in mean action potential duration
d3=(dmean3*mean_wt)/sd_wt
power10=tidy(power.t.test(n=NULL,sig.level = .05,power=0.8,sd=1,delta=d3))
#concatenate results
mean_difference=c(40,20,10)
cohen_d=rbind(d1,d2,d3)
sample_size=rbind(power40$n,power20$n,power10$n)
alpha=rbind(power40$sig.level,power20$sig.level,power10$sig.level)
power=rbind(power40$power,power20$power,power10$power)
lit_values=cbind(mean_difference,cohen_d,sample_size,alpha,power)
colnames(lit_values)=c("mean_difference","cohen_d","sample_size","alpha","power")
Q4.mean.difference.result.csv=write.csv(lit_values, "~/Desktop/Q4.mean.difference.result.csv")
print(lit_values)
## mean_difference cohen_d sample_size alpha power
## d1 40 2.0736445 4.828827 0.05 0.8
## d2 20 1.0368222 15.623711 0.05 0.8
## d3 10 0.5184111 59.386047 0.05 0.8
#experiment 1 40% change in mean action potential duration n=5
d1=(dmean1*mean_wt)/sd_wt
power40=tidy(power.t.test(n=5,sig.level = .05,power=NULL,sd=1,delta=d1))
#experiment 1 20% change in mean action potential duration n=5
d2=(dmean2*mean_wt)/sd_wt
power20=tidy(power.t.test(n=5,sig.level = .05,power=NULL,sd=1,delta=d2))
#experiment 1 10% change in mean action potential duration n=5
d3=(dmean3*mean_wt)/sd_wt
power10=tidy(power.t.test(n=5,sig.level = .05,power=NULL,sd=1,delta=d3))
#concatenate results
mean_difference=c(40,20,10)
cohen_d=rbind(d1,d2,d3)
sample_size=rbind(power40$n,power20$n,power10$n)
alpha=rbind(power40$sig.level,power20$sig.level,power10$sig.level)
power=rbind(power40$power,power20$power,power10$power)
lit_values=cbind(mean_difference,cohen_d,sample_size,alpha,power)
colnames(lit_values)=c("mean_difference","cohen_d","sample_size","alpha","power")
Q4.exp.1.result.csv=write.csv(lit_values, "~/Desktop/Q4.exp.1.result.csv")
print(lit_values)
## mean_difference cohen_d sample_size alpha power
## d1 40 2.0736445 5 0.05 0.8182358
## d2 20 1.0368222 5 0.05 0.3034952
## d3 10 0.5184111 5 0.05 0.1085460
#75% change in mean action potential
d_thesis=(thesis*mean_wt)/sd_wt
power75=tidy(power.t.test(n=NULL,sig.level = .05,power=0.8,sd=1,delta=d_thesis))
Q4.thesis.csv=write.csv(power75, "~/Desktop/Q4.thesis.csv")
print(power75)
## # A tibble: 1 x 5
## n delta sd sig.level power
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2.46 3.89 1 0.05 0.8