#=================== Apartado a)  ===================


library(PASWR)
attach(biomass)   
par(mfrow=c(1,2))
plot(Dn, PSA) # a clear non-linear relationship
abline(lsfit(Dn,PSA), col=2, lwd=2)
plot(log(Dn), log(PSA))   #the relation seems to be linear 
abline(lsfit(log(Dn),log(PSA) ), col=2, lwd=2)



model.DN<-lm(log(PSA)~log(Dn))
summary(model.DN)

par(mfrow=c(2,2))
plot(model.DN) 
 


#=================== Apartado b)  ===================

model.DNH<-lm(log(PSA)~log(Dn)+ H )
summary(model.DNH)

par(mfrow=c(2,2))
plot(model.DN) 
 

## The R^2 coefficient increases. (But the difference is negligible)
 

#=================== Apartado c) ===================

###### 1

summary(model.DNH)

###### 2

sigma2<-summary(model.DNH)$sigma^2
mcov<-summary(model.DNH)$cov.unscaled*sigma2
mcov

vcov(model.DNH)

###### 3

confint(model.DNH)

######  4

rsquared<-summary(model.DNH)$r.squared
rsquared.ad<-summary(model.DNH)$adj.r.squared

sigma2<-summary(model.DNH)$sigma^2
#sigma2<-sum(model.DNH$res^2)/39


######  5

win.graph()
par(mfrow=c(2,2),pty="s")
plot(model.DNH)

######  6

library(lmtest)
bptest(model.DNH)
#        studentized Breusch-Pagan test

#data:  model.DNH 
#BP = 7.626, df = 2, p-value = 0.02208




#######  7

shapiro.test(model.DNH$res)

#        Shapiro-Wilk normality test

#data:  model.DNH$res 
#W = 0.9563, p-value = 0.1088

shapiro.test(rstandard(model.DNH))
#        Shapiro-Wilk normality test
#
#data:  rstandard(model.DNH) 
#W = 0.9569, p-value = 0.1146

qqnorm(rstandard(model.DNH))
qqline(rstandard(model.DNH))

#######  8

a=model.DNH

win.graph()
plot(rstudent(a),type="n",xlab="",ylab="r_i^*")
text(rstudent(a))
abline(h=qt(0.025, a$df.residual-1))
abline(h=qt(0.975,a$df.residual-1))
title("c) Studentized Residuals")


# la 10, 11, 13, 15 




#######  9

lm.influence(model.DNH)
influence.measures(model.DNH)

summary(influence.measures(model.DNH))
 

### Otra forma

par(mfrow=c(2,3))
plot(model.DNH,which=1:6)


###### Mas manual 


iden<-function(y, a = 3, c=0.05)
{
        n <- length(y)
        oy <- order(abs(y))
        b<-y*c
        which <- oy[(n - a + 1):n]
        text(seq(1:n)[which], y[which]+b[which], as.character(which))
list(y=y,b=b)}

 # Cook distance and Dffits 


par(mfrow=c(1,2),pty="s")
cd.F<-cooks.distance(a)
plot(cd.F, ylab="Cooks Distance", ylim=c(0,0.8))
iden(cd.F, a=3)
crit.value<-qf(0.5, ncol(X), nrow(X)-ncol(X))
abline(h=crit.value, lty=2)

dffits.modelF<-dffits(a)
plot(dffits.modelF, ylab="Dffits", ylim=c(-1,1))
iden(dffits.modelF, a=3)
crit.value<-2*sqrt(ncol(X)/nrow(X))
abline(h=c(-crit.value, crit.value), lty=2)


#DFbetas Model F


par(mfrow=c(2,2))
dfbetas.modelF<-dfbetas(a)

plot(dfbetas.modelF[,1], ylab="dfbetas[,2]", ylim=c(-1,1))
iden(dfbetas.modelF[,1], a=3)
crit.value<-2/sqrt(nrow(X))
abline(h=c(-crit.value, crit.value), lty=2)


plot(dfbetas.modelF[,2], ylab="dfbetas[,2]", ylim=c(-1,1))
iden(dfbetas.modelF[,2], a=3)
crit.value<-2/sqrt(nrow(X))
abline(h=c(-crit.value, crit.value), lty=2)


plot(dfbetas.modelF[,3], ylab="dfbetas[,3]", ylim=c(-1,1))
iden(dfbetas.modelF[,3], a=3)
crit.value<-2/sqrt(nrow(X))
abline(h=c(-crit.value, crit.value), lty=2)

  

#=================== Apartado d) ===================

Dn<-seq(12.5,42.5,5)
H<-seq(10,40,5)
newdata<-data.frame(Dn,H)
predictions<-exp(predict.lm(model.DNH,newdata)+summary(model.DNH)$sigma^2/2)
predictions
                         
#=================== Apartado e)   ===================
                      
model.DNH2<-lm(PST~Dn+H,data=biomass) 
        

#=================== Apartado f)   ===================     

par(mfrow=c(2,2),pty="s")
plot(model.DNH2)
           
## There is a quadratic trend in the residuals   
      
#=================== Apartado g)   ===================     

model.DNH3<-lm(PST~Dn+I(Dn^2)+H,data=biomass) 
        
#=================== Apartado g)   ===================   
 
win.graph()                                                        
par(mfrow=c(2,2),pty="s")                                
plot(model.DNH3) 

## The quadratic trend has been mostly dissapeared, but there might be outliers.