#sample size
N=500

#Storage for the variables we will make
Av_qual_papers=rep(NA,N)
EV_Av_qual_papers=rep(NA,N)

#simulating the predictor variables
pubpress=rbinom(N,size=1,prob=0.6)
motiv = rbinom(N,size=1,prob=0.5) 

#junk for plotting colors later
colfilter = rep(NA,N)
colfilter[motiv==0]='blue'
colfilter[motiv==1]='deepskyblue'
colors1 <- colfilter

colfilter = rep(NA,N)
colfilter[pubpress==0]='pink'
colfilter[pubpress==1]='magenta'
colors2 <- colfilter

#sticking everything we have so far in a dataset
dataf=data.frame(Av_qual_papers,motiv,colors1,colors2,pubpress)

##simulating the dependent variable and sticking it in the dataset:
int=80 
B1=-8
B2=-30
B3=-20
error=rnorm(N,0,5)    

dataf[,"Av_qual_papers"] =int + B1*pubpress + B2*motiv + B3*pubpress*motiv + error

#calculate the expected values and also put this in our dataset
dataf[,"EV_Av_qual_papers"] =int + B1*pubpress + B2*motiv + B3*pubpress*motiv


##data analysis##
#maineffects=lm(Av_qual_papers ~ pubpress + motiv,data=dataf)
twowayinteractions=lm(Av_qual_papers ~  pubpress + motiv + pubpress*motiv,data=dataf)
summary(twowayinteractions)

#plots#
library('ggplot2')

plotje=ggplot(dataf, aes(x=pubpress, y=EV_Av_qual_papers, group=colors1, color=colors1))+ geom_point(size=2)+ geom_line(lwd=1.5) 
plotje+   labs(x="Publication Pressure", y="Average Quality Papers", 
               col="Group") + scale_color_manual(labels = c("Truthfinding", "Status"), values=c("blue","deepskyblue"))


plotje_rev=ggplot(dataf, aes(x=motiv, y=EV_Av_qual_papers, group=colors2, color=colors2))+ geom_point(size=2)+ geom_line(lwd=1.5) 
plotje_rev+   labs(x="Motivation", y="Average Quality Papers", 
               col="Group") + scale_color_manual(labels = c("high pressure", "low pressure"), values=c("pink","magenta"))


plot=ggplot(dataf, aes(x=pubpress, y=Av_qual_papers, group=colors1, color=colors1))+ geom_point(size=2)+ geom_line(lwd=1.5) 
plot+   labs(x="Publication Pressure", y="Average Quality Papers", 
               col="Group") + scale_color_manual(labels = c("Truthfinding", "Status"), values=c("blue","deepskyblue"))