rm(list=ls())

## Diamonds Price
## https://www.kaggle.com/datasets/shivam2503/diamonds

set.seed(87722)

D = read.csv("diamonds.csv",head=TRUE)
n = nrow(D)
I = sample(1:n,size=100,replace=FALSE)

D = D[I,]
D = D[,c(2,6,7,8)]
attach(D)
## pairs(D)
## names(D)

### simple Regression

pdf("Diamonds1.pdf")
## par(mfrow=c(2,2))
out = lm(price ~ carat)
plot(carat,price,xlab="",ylab="", pch=21,cex=2.,bg="pink",cex.axis=1.4)
abline(out,col="purple",lwd=4)
## print(summary(out))
title(xlab="Carats", line=2.5, cex.lab=1.8)
title(ylab="Price",  line=2.3, cex.lab=1.8)

### Loess
LO = loess(price ~ carat)
o=order(carat)
lines(carat[o],LO$fitted[o],col="red",lwd=3)
dev.off()

r = rstudent(out)
C = cooks.distance(out)

## Full regression

out1 = lm(price ~ carat + table + depth)
## summary(out1)

## Call:
## lm(formula = price ~ carat + table + depth)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2853.7  -870.8    64.7   705.0  5082.6 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11845.52    9753.76   1.214   0.2276    
## carat        8598.90     297.16  28.937   <2e-16 ***
## table         -35.97      75.95  -0.474   0.6369    
## depth        -204.83     116.99  -1.751   0.0832 .  
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 1424 on 96 degrees of freedom
## Multiple R-squared:  0.9045,	Adjusted R-squared:  0.9016 
## F-statistic: 303.2 on 3 and 96 DF,  p-value: < 2.2e-16


confint(out1)

##                2.5 %      97.5 %
## (Intercept) -7515.5545 31206.58593
## carat        8009.0490  9188.74974
## table        -186.7314   114.79923
## depth        -437.0465    27.38275

## Back to simple linear regression

pdf("ResidDiamonds.pdf")
plot(carat,rstudent(out),xlab="",ylab="", xaxt="n",pch=21,cex=2.2,
        bg="pink",cex.axis=1.5)
abline(h=0,col="red",lwd=4)
title(xlab="Carats", line=2, cex.lab=2)
title(ylab="Residuals", line=2.2, cex.lab=2)
dev.off()


pdf("CooksDiamonds.pdf")
C = cooks.distance(out)
r = rstudent(out)
plot(abs(r),C,xlab="",ylab="",xaxt="n",pch=21,cex=2.2,
     bg="pink",cex.axis=1.5)
title(xlab="Absolute Residuals",line=2, cex.lab=2)
title(ylab="Cook's Distances", line=2.2, cex.lab=2)
dev.off()


### Log transform both variables
pdf("RegrLogDiam.pdf")
logcarat = log(carat)
logprice = log(price)
plot(logcarat,logprice,xlab="",ylab="", pch=21,cex=2.2,
       bg="pink",cex.axis=1.4)
out1 = lm(logprice ~ logcarat)
abline(out1,col="purple",lwd=18)
## print(summary(out1))

LO2 = loess(logprice ~ logcarat)
o=order(logcarat)
lines(logcarat[o],LO2$fitted[o],col="red",lwd=3)
title(xlab="log(Carats)", line=2.5, cex.lab=2)
title(ylab="log(Price)",  line=2, cex.lab=2)
dev.off()

rlog = rstudent(out1)
Clog = cooks.distance(out1)

pdf("ResidLogDiam.pdf")
plot(logcarat,rlog,xlab="",ylab="",xaxt="n", pch=21,cex=2.2,
     bg="pink",cex.axis=1.5)
abline(h=0,col="red",lwd=4)
title(xlab="log(Carats)", line=2, cex.lab=2.2)
title(ylab="Residuals", line=2.2, cex.lab=2.2)
dev.off()

pdf("CooksLogDia.pdf")
plot(abs(rlog),Clog,xlab="",ylab="",xaxt="n",pch=21,cex=2,bg="pink")
title(xlab="Absolute Residuals",line=2, cex.lab=2)
title(ylab="Cook's Distances", line=2.2, cex.lab=2)
dev.off()


### PREDICTION using simple linear log transform. model
###            (out1=lm(logprice ~ logcarat))

### Interval predicting

## Confidence band for mu
pdf("ConfidenceMu.pdf")
range(logcarat)
## [1] -1.4696760  0.8197798
LC =seq(-1.4,0.8,length=100)
newD = data.frame(logcarat=seq(-1.4,0.8,length=100))
tmp = predict(out1, newdata=newD,interval="confidence",level=.95)
plot(LC,tmp[,1],lwd=3,ylim=c(min(tmp),max(tmp)),type="l",xlab="",
     ylab="",col="purple",cex.axis=1.5)
polygon(c(rev(LC),LC),c(rev(tmp[,3]),tmp[,2]),density=NULL,angle=90,
border=NA,col="thistle",lty=1) 
lines(LC,tmp[,1],col="purple",lty=1,lwd=1)
lines(LC,tmp[,2],col="purple",lty="dashed",lwd=2)
lines(LC,tmp[,3],col="purple",lty="dashed",lwd=2)
title(xlab="log(Carats)", line=2.5, cex.lab=2)
title(ylab=expression(mu(x)),  line=2, cex.lab=2)
dev.off()

### Prediction bands for \hat(Y) 

pdf("PredictYDiam.pdf")
tmp = predict(out1, newdata=newD,interval="prediction",level=.95)
plot(LC,tmp[,1],lwd=3,ylim=c(min(tmp),max(tmp)),type="l",col="purple",
     xlab="",ylab="",cex.axis=1.3)
polygon(c(rev(LC),LC),c(rev(tmp[,3]),tmp[,2]),density=NULL,angle=90,
        border=NA,col="thistle",lty=1)      
lines(LC,tmp[,1],col="purple",lty=4,lwd=3)
lines(LC,tmp[,2],col="purple",lwd=1)
lines(LC,tmp[,3],col="purple",lwd=1)
title(xlab="log(Carats)", line=2.7, cex.lab=2.1)
title(ylab=expression(hat(Y)),  line=1.9, cex.lab=1.8)
dev.off()



## plot(table,r,xlab="",ylab="", pch=21,cex=1.1,bg="pink")
## title(xlab="Table", line=2.5, cex.lab=1.2)
## title(ylab="Residuals multiple", line=2.5, cex.lab=1.2)