##### Demonstration on reading regression outputs and residuals. F. Borghesi 27 Oct 2021


### 1 CATEGORICAL VARIABLE (ANOVA)
n<-200   #set sampel size
cat1<-rnorm(n,10,3) # creates data for category 1 with mean = 10 and sd= 3
cat2<-rnorm(n,8,3) # creates data for category 1 with mean = 8 and sd= 3
cat3<-rnorm(n,7.5,3)  # creates data for category 1 with mean = 7.5 and sd= 3

data<-data.frame(VarX=c(rep("Cat1",n),rep("Cat2",n),rep("Cat3",n)),
             VarY=c(cat1,cat2,cat3))  #just joins the above on a data frame

model<-lm(VarY~VarX,data=data)  # creates lm model (with single categorical variable)
summary(model)
plot(data$VarX,data$VarY)

hist(data$VarY[data$VarX=="Cat1"])
hist(residuals(model)[data$VarX=="Cat1"])

shapiro.test(data$VarY[data$VarX=="Cat1"])   # checks normality for cat1 data
shapiro.test(residuals(model)[data$VarX=="Cat1"])   # checks normality for cat1 residuals

### 1 QUANTITATIVE VARIABLE (SIMPLE REGRESSION)
n<-400
int<-1.4  # sets an intercept
slp<-1.21
x<-runif(n,0,23) 
y<-int+slp*x+rnorm(n,0,4) # sets linear response with normal deviations!!!
plot(x,y)


model<-lm(y~x)
summary(model)
abline(model)

hist(y)
hist(residuals(model))
shapiro.test(y)
shapiro.test(residuals(model))


### 1 QUANTITATIVE VARIABLE AND 1 CATEGORICAL VARIABLE (ANCOVA)
n<-400
int1<-1.4
int2<-2.8
slp1<-1.21
slp2<-2.16

x1<-runif(n,0,23)
x2<-runif(n,0,23)
y1<-int1+slp1*x1+rnorm(n,0,4)
y2<-int2+slp2*x2+rnorm(n,0,4)

data<-data.frame(VarC=c(rep("Cat1",n),rep("Cat2",n)),VarX=c(x1,x2),
                 VarY=c(y1,y2))

plot(data$VarX,data$VarY,col=factor(data$VarC))

model<-lm(VarY~VarC+VarX,data=data)
summary(model)

plot(data$VarX,data$VarY,col=factor(data$VarC))
abline(model$coefficients[1],model$coefficients[3])
abline(model$coefficients[1]+model$coefficients[2],model$coefficients[3],col="red")


model2<-lm(VarY~VarC*VarX,data=data)
summary(model2)

plot(data$VarX,data$VarY,col=factor(data$VarC))
abline(model2$coefficients[1],model2$coefficients[3])
abline(model2$coefficients[1]+model2$coefficients[2],
       model2$coefficients[3]+model2$coefficients[4],col="red")

hist(data$VarY[data$VarC=="Cat1"])
hist(residuals(model2)[data$VarC=="Cat1"])

shapiro.test(data$VarY[data$VarC=="Cat1"])
shapiro.test(residuals(model2)[data$VarC=="Cat1"])

car::leveneTest(data$VarY~data$VarC)
car::leveneTest(residuals(model2)~data$VarC)