# R script for multiple regression using lm() function.
# Data in lungpressure.txt:  Y is an invasive method for
# determining arterial lung pressure (potentially risky to
# patient). Predictor variables: X1=emptying rate of blood
# into pumping chamber of heart (radionuclide imaging),
#  X2=ejection rate of blood pumped out of heart to lungs
# (radionuclide imaging), X3=a blood gas concentration. 
#
# Data from A. T. Marmor et al. (1986) Chest 89:64-69
# as adapted by Kutner et al. (2004) Applied linear
# regression models, fourth edition . McGraw-Hill.
#
# All subsets regression required installation of the "leaps"
# package.

data=read.table("lungpressure.txt",header=TRUE)
attach(data)
head(data)

# Calculate & print regression model 1.
results=lm(y~x1+x2+x3,data)
summary(results)
residuals(results)
fitted(results)
AIC(results)

# Other useful functions 
coefficients(results) # model coefficients
confint(results, level=0.95) # CIs for model parameters 
fitted(results) # predicted values
anova(results) # anova table 
vcov(results) # covariance matrix for model parameters 
influence(results) # regression diagnostics 


layout(matrix(1:4,2,2))
plot(results)                       #  Diagnostic plots for evaluating
                                    #  the regression model.

# Calculate & print regression model 1 without intercept.
results=lm(y~0+x1+x2+x3,data)
summary(results)

# All subsets regression.
library(leaps)
all.mods=regsubsets(y~x1+x2+x3+x1*x2+x1*x3+x2*x3+x1*x2*x3,
           data)

# view results 
summary(all.mods)
# plot a table of models showing variables in each model.
# models are ordered by the selection statistic.
plot(all.mods,scale="Cp")  # Cp orders models identical to AIC.

# Calculate & print regression model 2:  x1, x2, x1*x2 
#  (x1, x2, plus interaction between x1 and x2).
results=lm(y~x1+x2+x1*x2,data)
summary(results)

layout(matrix(1:4,2,2))
plot(results)                       #  Diagnostic plots for evaluating
                                    #  the regression model. 

detach(data)