# R code for chapter 59 Multiple Linear regression chapter
# Book details:
# Health Science Statistics using R and R Commander
# by Robin Beaumont
# Paperback - ISBN 9781907904318
# http://www.amazon.co.uk/Health-Science-Statistics-using-Commander/dp/190790431X
# or direct from publisher at: 
# http://www.scionpublishing.com/
# book website: http://www.robin-beaumont.co.uk/rbook
###########################################################
# Multiple Linear regression chapter

# get file locally ignore:
library(foreign)
# campbell_mult_reg <- 
#   read.spss("D:/web_sites_mine/HIcourseweb new/book2data/campbell_mlr_p12.sav",
#    use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)
# names(campbell_mult_reg)
###### get file from web

campbell_mult_reg <- 
  read.spss("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/campbell_mlr_p12.sav",
   use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)
names(campbell_mult_reg)

# code generated by R Commander when selecting correlation matrix, partial correlations:
partial.cor(campbell_mult_reg[,c("age","deadsp","height")], 
tests=TRUE,   use="complete")
 
 ####
install.packages("ppcor", dependencies=TRUE)
library(ppcor)
# partial correlations using pcor.test() 

# part correlations sp.test()
pcor.test(campbell_mult_reg$deadsp, campbell_mult_reg$age, campbell_mult_reg$height, method="kendall")
pcor.test(campbell_mult_reg$deadsp, campbell_mult_reg$age, campbell_mult_reg$height, method="spearman")
partial_Rsqr <- pcor.test(campbell_mult_reg$deadsp, campbell_mult_reg$age, campbell_mult_reg$height)
partial_Rsqr
# [row, column]

spcor(campbell_mult_reg[,c(2,3,5)])
pcor(campbell_mult_reg[,c(2,3,5)])
########  equivalent to page 13 in chapter
reg_height_age <- lm(deadsp~age+height, data=campbell_mult_reg)
 summary(reg_height_age)
########  equivalent to page 14 in chapter
 #  for standardised betas edited above to: 
reg_height_age2 <- lm( scale(deadsp)~ scale(age) + scale(height), data=campbell_mult_reg)
 summary(reg_height_age2)
########  equivalent to page 14 in chapter
 #  VIF's 
 vif(reg_height_age)
 sqrt(3.73.8)
 ########  equivalent to page 15 in chapter
 #  confidence intervals 
 library(MASS, pos=4)
Confint(reg_height_age, level=0.95)

 ########  equivalent to page 15 in chapter
 #  anova tables
anova(reg_height_age)
Anova(reg_height_age, type="I") # produces error must be II or III
Anova(reg_height_age, type="II")
Anova(reg_height_age, type="III")

###### create height only model
########  equivalent to page 16 in chapter
height_only <- lm(deadsp~height, data=campbell_mult_reg)
summary(height_only)

######## comparing models
anova( reg_height_age, height_only)
######
###### AIC and BIC
########  equivalent to page 17 in chapter
AIC(height_only)
BIC(height_only)
AIC(reg_height_age)
BIC(reg_height_age)
#################### next section in chapter
##########  Two predictors, one continuous one nominal
## model1 <- deadsp~ asthma + height
model1 <- lm(deadsp ~ asthma +height, data=campbell_mult_reg)
summary(model1)
 ########  equivalent to page 20 in chapter
Draw the scatterplot - generated by R commander:
scatterplot(deadsp~height | asthma, reg.line=FALSE, smooth=FALSE, 
  spread=FALSE, id.method='mahal', id.n = 2, boxplots=FALSE, span=0.5, 
  by.groups=TRUE, data=campbell_mult_reg)
  # add the lines
  abline(a=-46.29, b=0.845, col = "blue", lty = 1)
abline(a=-63.1, b=0.845, col = "red", lty = 3)
# and get a confidence interval:
Confint(model1, level=0.95)

#now get the anova tables
anova(model1) # give total for asthma first
Anova(model1, type="II")
Anova(model1, type="III")
 ########  equivalent to page 22 in chapter
 # now non parallel lines - code generated by R commander:
 
 scatterplot(deadsp~height | asthma, reg.line=lm, smooth=FALSE, spread=FALSE,
   id.method='mahal', id.n = 2, boxplots=FALSE, span=0.5, by.groups=TRUE, 
  data=campbell_mult_reg)
 ######### 
  ########  equivalent to page 23 in chapter
 
# ANCOVA + interaction to assess it
########  equivalent to page 23 chapter
model2<- lm(deadsp ~ asthma * height, data=campbell_mult_reg)
summary(model2)
# basic diagnostic plots
oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
plot(model2)
par(oldpar)


anova(model2) # give total for asthma first
Anova(model2, type="II")
Anova(model2, type="III")
model2_temp<- lm(deadsp ~ height * asthma  , data=campbell_mult_reg)
summary(model2_temp)
anova(model2_temp) # give total for height first
Anova(model2_temp, type="II")
Anova(model2_temp, type="III")

###  equivalent to page 24 chapter
#### centering for vif's
model2_st<- lm(scale(deadsp) ~ asthma * scale(height), data=campbell_mult_reg)
summary(model2_st)
model2_cent<- lm((deadsp - mean(deadsp)) ~ asthma * scale(height-mean(height)), data=campbell_mult_reg)
summary(model2_cent)

##############
# relative importance of inputs
########  equivalent to page 24 in chapter
install.packages("relaimpo", dependencies=TRUE)
library(relaimpo)
average_partial_estimates <- calc.relimp(model2, type = "lmg")
average_partial_estimates

average_partial_estimates2 <- calc.relimp(model2_temp, type = "lmg",rela=FALSE)
average_partial_estimates2

#################
# to produce the tableplot for collinearity metrics
# first need to perturb package
#‘‘colldiag’’ calculates condition numbers and
#variance decomposition proportions to test for collinearity and uncover its sources
#colldiag(mod, scale = TRUE, center = FALSE, add.intercept = TRUE)
# Function 'colldiag' with default arguments (standardizes the data).
#coll.1 <- colldiag(mod = mod.1, scale = TRUE, center = FALSE, add.intercept = TRUE)
# coll.1
# Function 'colldiag' with optional values for arguments (unstandardized).
# coll.2 <- colldiag(mod = mod.1, scale = FALSE, center = FALSE, add.intercept = TRUE)
# coll.2
install.packages("perturb")
library(perturb)
# THE coldiag() function only works with numeric inputs
# the first model below produces a error: Error in v^2 : non-numeric argument to binary operator
# in other words it can't cope with factors
collinearity_diags1 <- colldiag(model2, scale = TRUE, center = FALSE, add.intercept = TRUE)
# this works
collinearity_diags1 <- colldiag(reg_height_age)
collinearity_diags1
collinearity_diags2 <- colldiag(mod = reg_height_age, scale = FALSE, center = FALSE, add.intercept = TRUE)
collinearity_diags2

########
install.packages("tableplot")
library(tableplot)

cd <- colldiag(model2, add.intercept=FALSE, center=TRUE)
print(cd, fuzz=.3)
tableplot(collinearity_diags1)

################## contrasts
reg_factor3levels <- lm(deadsp~age+height +  condition, data=campbell_mult_reg)
 summary(reg_factor3levels)
contrasts(campbell_mult_reg$condition)
# gives: asthma=0,0; bronchitis=1,0; normal=0,1
contrasts (campbell_mult_reg$condition <- contr.treatment(3, base = 3)
reg_factor3levelsb <- lm(deadsp~age+height +  condition, data=campbell_mult_reg)
 summary(reg_factor3levelsb)

colnames(contrasts(campbell_mult_reg$condition))<- c("asthma_v_normal", "Bronchitis_v_normal")
contrasts(campbell_mult_reg$condition)
reg_factor3levelsb <- lm(deadsp~age+height +  condition, data=campbell_mult_reg)
 summary(reg_factor3levelsb)
 scale
 reg_factor3levelsb <- lm(deadsp~  scale(age) +  scale(height) +  condition, data=campbell_mult_reg)
 summary(reg_factor3levelsb)

tapply(campbell_mult_reg$deadsp, campbell_mult_reg$condition, mean)

reg_factor3levelsc <- lm(deadsp~ condition, data=campbell_mult_reg)
 summary(reg_factor3levelsc)
 
 ################ predcting new values
model2<- lm(deadsp ~ asthma * height, data=campbell_mult_reg)
 predict(model2, data.frame(asthma="yes", height=120))
 # and with prediction or conf interval
  predict(model2, data.frame(asthma="yes", height=120), interval="pred")
   predict(model2, data.frame(asthma="yes", height=120), interval="conf", level=.95) 
   #################### additional diagnostics
   library(car)
   residualPlots(model2)
#############################################
# basic diagnostic plots
# need 4 plots on a single page
oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
plot(reg_height_age)
par(oldpar)
influence.measures(reg_height_age)
###
Influence measures of
#         lm(formula = deadsp ~ age + height, data = campbell_mult_reg) :
#
#     dfb.1_  dfb.age dfb.hght   dffit cov.r   cook.d    hat
#1   0.59513  0.22331 -0.47676  0.6659 1.538 0.149159 0.3320
#2  -0.25431 -0.03930  0.17605 -0.3311 1.516 0.038541 0.2256
#3  -0.10209 -0.02961  0.07556 -0.1344 1.496 0.006508 0.1550
#4  -0.23107 -0.13900  0.20243 -0.2877 1.340 0.028843 0.1474
#5   0.02429  0.01315 -0.02039  0.0328 1.475 0.000390 0.1215
#6  -0.11334 -0.51645  0.33711  0.7978 0.554 0.166451 0.1290
#7   0.03533  0.07572 -0.06044 -0.0972 1.550 0.003423 0.1728
#8   0.10870  0.10619 -0.14169 -0.3697 0.938 0.043397 0.0787
#9   0.21538  0.20373 -0.25635 -0.4373 0.945 0.060345 0.1023
#10 -0.08818  0.02190  0.07076  0.3565 1.012 0.041254 0.0876
#11 -0.20590 -0.24144  0.25227  0.2951 1.687 0.031025 0.2765
#12  0.40800  0.37914 -0.46084 -0.5749 1.108 0.106377 0.1878
#13 -0.00917  0.00144  0.00652  0.0209 1.507 0.000158 0.1390
#14 -0.25166  0.28589  0.07985  0.8518 0.773 0.206700 0.1926
#15 -0.14989 -0.78909  0.41966 -0.9822 3.255 0.335130 0.6522

Also 
paste(resid(reg_height_age), rstandard(reg_height_age),rstudent(reg_height_age))

result_table <- cbind( 'raw res' = resid(reg_height_age), 
'standardised res' = rstandard(reg_height_age), 
'studentized res' = rstudent(reg_height_age) ) 
result_table

######## PRESS statistic
# calculating PRESS statistic (for model validation)
# Computation of Allen’s PRESS statistic for an lm object: PRESS(lm) 

# from http://www.stat.sc.edu/~hitchcock/diagnosticRexample704.txt 
PRESS_statistic <- sum( (resid(reg_height_age)/(1-hatvalues(reg_height_age)))^2 )
print(paste("PRESS statistic= ", PRESS_statistic))
# [1] 3050.161
# alternatives PRESS in MPV package

install.packages("MPV", dependencies=TRUE)
library(MPV)
PRESS(reg_height_age)
# [1] 3050.161
detach("package:MPV", unload=TRUE)

### alternative qpcR PRESS() - works for glm as well as lm()
# produces the residuals the sum statistic and the P-square value (i.e. aka r squared?)
library(qpcR)
PRESS(reg_height_age)
# or
qpcR::PRESS(reg_height_age)
######################################
## appendix mahalanobis distances
## many different packages provide this - several varieties
## using 
## Robust Mahalanobis distance. 
# here you use the original data not the lm() object
# Using library(chemometrics) and the  Moutlier() function; 
# note: several dependent packages
install.packages("chemometrics", dependencies=TRUE)
# Get the robust Mahalanobis
# distances for each case of our current data; the plot contains both robust (rd) 
# and classical (md) Mahalanobis distances.
# use the campbell_mult_reg  dataset  but need to only include appropriate columns
# names(campbell_mult_reg)
# [1] "no"         "deadsp"     "height"     "asthma"     "age"       
# [6] "bronchitis" "normal"     "ht_ath"     "condition" 
# reminder: 
# reg_height_age <- lm(deadsp~age+height,
# model1 <- lm(deadsp ~ asthma +height,
# model2<- lm(deadsp ~ asthma * height,

attach(campbell_mult_reg)
subset_height_age <- data.frame(deadsp, age, height)
subset_height_age
subset_height_age_only <- data.frame(age, height)
subset_height_age_only
detach(campbell_mult_reg)
# rm(subset.1)
attach(subset_height_age)
library(chemometrics)
trad_m2 <- Moutlier(subset_height_age_only,quantile=0.975,plot=FALSE)$md
trad_m2
summary(trad_m2)
# $rd = robust 
robust_m2 <- Moutlier(subset_height_age_only,quantile=0.975,plot=TRUE)$rd
robust_m2
summary(robust_m2)
# Q to Q plot of the robust Mahalanobis distance (used to look for outliers). 
# qqPlot() in the car package
library(car)
par(mfrow=c(1,1))
qqPlot(robust_m2, dist= "norm", line="robust")
legend("topleft", "Robust Mahalanobis Distances", pch=1, title="Quantile-Quantile Plot", col="red")


# If desired; you can get the classic Mahalanobis distances for each case 
# of our current data.
# $md = robust 
trad_m <- Moutlier(subset_height_age,quantile=0.975,plot=FALSE)$md
summary(trad_m)
# qqPlot() in the car package
par(mfrow=c(1,1))
qqPlot(trad_m, dist= "norm", line="robust")
legend("topleft", "Classic Mahalanobis Distances", pch=1, title="Quantile-Quantile Plot", col="red")
par(oldpar)

## Comparison of Classic & Robust Mahalanobis Distances.

scatterplot(robust_m, trad_m, smooth=TRUE, ellipse=TRUE, span =1,
     reg.line=lm, boxplots='xy')

# Cases with extreme Mahalanobis distances (upper right) are likely to be 
# true multivariate outliers.
################ first principles from:
# https://stat.ethz.ch/R-manual/R-devel/library/stats/html/mahalanobis.html
# adapted
Sx <- cov(subset_height_age)
obs <- nrow(subset_height_age)
pedictors <- 2
D2 <- mahalanobis(subset_height_age, colMeans(subset_height_age), Sx)
plot(density(D2, bw = 0.5),
     main=paste("Squared Mahalanobis distances, n=", obs,"  p=2") )
	 rug(D2)
	 
qqplot(qchisq(ppoints(obs), df = 2), D2,
       main = expression("Q-Q plot of Mahalanobis" * ~D^2 *
                         " vs. quantiles of" * ~ chi[3]^2))
abline(0, 1, col = 'gray')