#No ethical control
N=500
Av_qual_papers=rep(NA,N)
EV_Av_qual_papers=rep(NA,N)
pubpressm=rep(0,N)
pubpressl=pubpressm

pubpressm[1:round(1/3*N)]= rep(1,round(1/3*N))     #0/1  geen druk medium druk
pubpressl[(round(1/3*N)+1):(2*(round(1/3*N)))]=rep(1,round(1/3*N))      #0/1  geen druk veel druk
motiv = rbinom(N,size=1,prob=0.5)   #0/1 waarheid vs status
ethictrl = rbinom(N,size=1,prob=0.5)    # 0/1 (geen controle veel controle)

pubpress=pubpressm
pubpress[pubpressl==1]=2

colfilter = rep(NA,N)
colfilter[motiv==0&ethictrl==1]='blue'
colfilter[motiv==1&ethictrl==1]='deepskyblue'
colfilter[motiv==0&ethictrl==0]='red'
colfilter[motiv==1&ethictrl==0]='orange'
colors <- colfilter
dataf=data.frame(Av_qual_papers,pubpressm,pubpressl,motiv,ethictrl,colors,pubpress)



int=80
B1=-8
B2=-50
B3=-40
B4=10
B5=-7       
B6=20
B7=25 
B8=5      
B9=40
B10=5
B11=-35

error=rnorm(N,0,5)       

dataf[,"EV_Av_qual_papers"] =int + B1*pubpressm +B2*pubpressl+ B3*motiv + B4*ethictrl+ B5*pubpressm*motiv+B6*pubpressl*motiv  + B7*ethictrl*motiv+ B8*ethictrl*pubpressm + B9*ethictrl*pubpressl+ B10*pubpressm*motiv*ethictrl + B11*pubpressl*ethictrl*motiv

dataf[,"Av_qual_papers"] =int + B1*pubpressm +B2*pubpressl+ B3*motiv + B4*ethictrl+ B5*pubpressm*motiv+B6*pubpressl*motiv  + B7*ethictrl*motiv+ B8*ethictrl*pubpressm + B9*ethictrl*pubpressl+ B10*pubpressm*motiv*ethictrl + B11*pubpressl*ethictrl*motiv + error

######
maineffects=lm(Av_qual_papers ~ pubpressm + pubpressl + motiv + ethictrl,data=dataf)
twowayinteractions=lm(Av_qual_papers ~ pubpressm + pubpressl + motiv + motiv*ethictrl +ethictrl+pubpressm*motiv+pubpressm*ethictrl+ethictrl+pubpressl*motiv+pubpressl*ethictrl,data=dataf)
threewayinteractions=lm(Av_qual_papers ~ pubpressm + pubpressl + motiv + motiv*ethictrl +ethictrl+pubpressm*motiv+pubpressm*ethictrl+ethictrl+pubpressl*motiv+pubpressl*ethictrl
 +pubpressm*ethictrl*motiv + pubpressl*ethictrl*motiv  ,data=dataf)

####################3DPlot############################
# library
library(rgl)

# This is to output a rgl plot in a rmarkdown document. Note that you must add webgl=TRUE, results='hide' in the chunck header
#library(knitr)
#knit_hooks$set(webgl = hook_webgl)

# Plot
par(mar=c(0,0,0,0))
plot3d( 
  y=dataf$pubpress, x=dataf$Av_qual_papers, z=dataf$ethictrl, 
  col = dataf$colors, 
  type = 's', 
  radius = 0.5,
  ylab="Publication Pressure", zlab="Ethical Control", xlab="Average Quality Papers",
  xlim=(c(0,100)))

#writeWebGL( filename="HtmlWidget/3dscatter.html" ,  width=600, height=600)

library('ggplot2')

ggplot(dataf, aes(x=pubpress, y=EV_Av_qual_papers, group=motiv,color=colors))+ geom_point() 

plotje=ggplot(dataf, aes(x=pubpress, y=EV_Av_qual_papers, group=colors, color=colors))+ 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("Ethical Control, Truthfinding", "Ethical Control, Status", "No Ethical Control, Status","No Ethical Control, Truthfinding"), values=c("blue","deepskyblue", "red", "orange"))

datafnoethc=dataf[-c(which(dataf$ethictrl==1)),]
plotje=ggplot(datafnoethc, aes(x=pubpress, y=EV_Av_qual_papers, group=colors, color=colors))+ 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("No Ethical Control, Status","No Ethical Control, Truthfinding"), values=c("red", "orange"))

datafethc=dataf[-c(which(dataf$ethictrl==0)),]
plotje=ggplot(datafethc, aes(x=pubpress, y=EV_Av_qual_papers, group=colors, color=colors))+ 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("Ethical Control, Truthfinding","Ethical Control, Status"), values=c("blue", "deepskyblue"))

 #Av_qual_papers = int + B1*pubpressm +B2*pubpressl+ B3*motiv + B4*ethictrl 
#+ B5*pubpressm*motiv+B6*pubpressl*motiv  + B7*ethictrl*motiv+ B8*ethictrl*pubpressm + B9*ethictrl*pubpressl
#+ B10*pubpressm*motiv*ethictrl + B11*pubpressl*ethictrl*motiv


####
#No pressure 00, truthfinding 0, no control 0.
#Av_qual_papers = int + B1*0 +B2*0+ B3*0 + B4*0 
#+ B5*0+ B6*0  + B7*0+ B8*0 + B9*0
#+ B10*0 + B11*0
#EV=int
####
####
#Medium pressure 10, truthfinding 0, no control 0.
#Av_qual_papers = int + B1*1 +B2*0+ B3*0 + B4*0 
#+ B5*1*0+B6*0*0  + B7*0*0+ B8*0*1 + B9*0*0
#+ B10*1*0*0 + B11*0*0*0
#EV=int+B1
####
####
#large pressure 01, truthfinding 0, no control 0.
#Av_qual_papers = int + B1*0 +B2*1+ B3*0 + B4*0 
#+ B5*0*0+B6*1*0  + B7*0*0+ B8*0*0 + B9*0*1
#+ B10*0*0*0 + B11*1*0*0
#EV=int+B2
####
####
#No pressure 00, status 1, no control 0.
#Av_qual_papers = int + B1*0 +B2*0+ B3*1 + B4*0 
#+ B5*0*1+ B6*0*1  + B7*0*1+ B8*0*0 + B9*0*0
#+ B10*0*1*0 + B11*0*0*1
#EV=int+B3
####
####
#Medium pressure 10, status 1, no control 0.
#Av_qual_papers = int + B1*1 +B2*0+ B3*1 + B4*0 
#+ B5*1*1+ B6*0*1 + B7*0*1+ B8*0*1 + B9*0*0
#+ B10*1*1*0 + B11*0*0*1
#EV=int+B1+B3+B5  
####
####
#large pressure 01, status 1, no control 0.
#Av_qual_papers = int + B1*0 +B2*1+ B3*1 + B4*0 
#+ B5*0*1+ B6*1*1  + B7*0*1+ B8*0*0 + B9*0*1
#+ B10*0*1*0 + B11*1*0*1
#EV=int+B2+B3+B6 


####
####
#No pressure 00, truthfinding 0, control 1.
#Av_qual_papers = int + B1*0 +B2*0+ B3*0 + B4*1 
#+ B5*0*0+B6*0*0  + B7*1*0+ B8*1*0 + B9*1*0
#+ B10*0*0*1 + B11*0*1*0
#EV=int+B4
####
####
#Medium pressure 10, truthfinding 0, control 1.
#Av_qual_papers = int + B1*1 +B2*0+ B3*0 + B4*1 
#+ B5*1*0+B6*0*0  + B7*1*0+ B8*1*1 + B9*1*0
#+ B10*1*0*1 + B11*0*1*0
#EV=int+B1+B4+B8
####
####
#large pressure 01, truthfinding 0,  control 1.
#Av_qual_papers = int + B1*0 +B2*1+ B3*0 + B4*1 
#+ B5*0*0+B6*1*0  + B7*1*0+ B8*1*0 + B9*1*1
#+ B10*0*0*1 + B11*1*1*0
#EV=int+B2+B4+B9
####
####
####
#No pressure 00, status 1, control 1.
#Av_qual_papers = int + B1*0 +B2*0+ B3*1 + B4*1 
#+ B5*0*1+B6*0*1  + B7*1*1+ B8*1*0 + B9*1*0
#+ B10*0*1*1 + B11*0*1*1
#EV=int+B3+B4+B7
####
####
#Medium pressure 10, status 1, control 1.
#Av_qual_papers = int + B1*1 +B2*0+ B3*1 + B4*1 
#+ B5*1*1+B6*0*1  + B7*1*1+ B8*1*1 + B9*1*0
#+ B10*1*1*1 + B11*0*1*1
#EV=int+B1+B3+B4+B5+B7+B8+B10
####
####
#large pressure 01, status 1,  control 1.
#Av_qual_papers = int + B1*0 +B2*1+ B3*1 + B4*1 
#+ B5*0*1+B6*1*1  + B7*1*1+ B8*1*0 + B9*1*1
#+ B10*0*1*1 + B11*1*1*1
#EV=int+B2+B3+B4+B6+B7+B9+B11
####