rm(list=ls())
setwd("/Users/eciliberto/Documents/corsi_statistica/R/R_definitivo/Introduction_R/data")  ###define your working directory (the folder where you saved data)
data<-read.table("fev.csv",head=T,sep="")
dim(data)
# We can get basic descriptives for the entire data set by using summary
names(data)
head(data)
summary(data)
str(data)

hist(data$fev,c(0,1,2,3,4,5,6),freq=F,main="Histogram",xlab="Fever",ylim=c(0,0.6))
m<-mean(data$fev)
std<-sqrt(var(data$fev))
x<-seq(min(data$fev), max(data$fev), length = 40)
y<-dnorm(x, mean=m, sd=std)
lines(x,y,col="red", lwd=2)

plot(data$age,data$fev,xlab="age",ylab="Fev")
abline(lm(data$fev~data$age))
mylinear.age<-lm(fev~age,data=data)
summary(mylinear.age)
confint(mylinear.age)
data$cage<-data$age-mean(data$age)
mylinear.age<-lm(fev~cage,data=data)
summary(mylinear.age)
confint(mylinear.age)

plot(data$height,data$fev,xlab="Height",ylab="Fev")
abline(lm(data$fev~data$height))
data$cheight<-data$height-mean(data$height)
mylinear.height<-lm(fev~cheight,data=data)
summary(mylinear.height)
confint(mylinear.height)

plot(data$sex,data$fev,xlab="Sex",ylab="Fev")
mylinear.sex<-lm(fev~sex,data=data)
summary(mylinear.sex)
confint(mylinear.sex)

mylinear.smoke<-lm(fev~smoke,data=data)
summary(mylinear.smoke)
confint(mylinear.smoke)

data$age.cat<-cut(data$age,breaks=c(0,10,15,20),right=F)
table(data$age.cat)
mylinear.agecat<-lm(fev~age.cat,data=data)
summary(mylinear.agecat)
confint(mylinear.agecat)

mylinear<-lm(fev~cage+cheight+sex+smoke,data=data)
summary(mylinear)
confint(mylinear)

qqnorm(residuals(mylinear))
shapiro.test(residuals(mylinear))
hist(residuals(mylinear))
plot(fitted(mylinear),residuals(mylinear),xlab="valori stimati",ylab="Residui")
abline(h=0)

newdata1<- with(data, data.frame(cage=0,cheight=0,sex=c("female","male","female","male"),smoke=c("current smoker","current smoker","non-current smoker","non-current smoker")))
newdata1
newdata2 <- cbind(newdata1, predict(mylinear, newdata = newdata1,se.fit = TRUE))
newdata2
pred.w.clim <- predict(mylinear, newdata1, interval = "confidence")
cbind(newdata2,pred.w.clim)