SPECIFIC SPLUS COMMANDS FOR BASIC STATISTICS AND PLOTS

postscript(file="plot.ps",horizontal=F)   ## Plots directed to file "plot.ps"
par(mfrow=c(3,1))

## Read Single Continuous Variable Data, obtain basic information

decaydat <- scan("decay-times.dat")   # read single variable data from a file
x1 <- c(429,430,430,431,436,437,440,441,445,446,447)   # enter single variable
x2 <- c(418,421,421,422,425,427,431,434,437,439)
summary(decaydat)       # provides basic summary statistics for one sample data
mean(decaydat)
var(decaydat)
hist(decaydat)          # provides histogram for one sample data
stem(decaydat)          # provides stem-and-leaf plot for one sample data
boxplot(decaydat)       # provides boxplot for one sample data
qqnorm(decaydat)        # provides normal probability plot for one sample data
boxplot(x1,x2)          # provides comparative boxplots for 2 samples
qqplot(x1,x2)           # provides quantile-quantile plot for 2 samples
u <- sqrt(decaydat)
qqnorm(u)
u <- log(decaydat)
qqnorm(u)
qqline(u)
hist(u)
## Construct a Q-Q Plot to Check Data Relative to Exponential Distribution
ps <- (c(1:n) - 0.5)/n          ## Values  (i - 0.5)/n
expscore <- - log(1-ps)         ## Values  F^{-1}[(i - 0.5)/n]
plot(expscore,sort(decaydat))


SPECIFIC SPLUS COMMANDS FOR LINEAR REGRESSION ANALYSIS AND ANOVA

postscript(file="plot.ps",horizontal=F)    ## Plots directed to file "plot.ps"
par(mfrow=c(3,1))

## Read in Data as Matrix, obtain basic information on pairwise relations

tobaccodat <- matrix(scan("/u/r/e/reinsel/stat333/tobacco.dat"),ncol=7,byrow=T)
dimnames(tobaccodat) <- list(NULL, c('x1','x2','x3','x4','x5','x6','y'))
pairs(tobaccodat)        # obtain scatterplot matrix of y and x variables
cor(tobaccodat)          # obtain correlation matrix of y and x variables
y <- tobaccodat[,7]
x1 <- tobaccodat[,1]
x2 <- tobaccodat[,2]
x3 <- tobaccodat[,3]
x4 <- tobaccodat[,4]
x5 <- tobaccodat[,5]
x6 <- tobaccodat[,6]

## Specify and Estimate Linear Model by LS

model1 <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6)  # estimate linear model by LS
summary(model1)                # display results from LS estimation of model
anova(model1)                  # display seq. SS for all terms in regression

model2 <- lm(y ~ x1 + x2 + x6)        # estimate 2nd linear (sub)model by LS
summary(model2)                # display results from LS estimation of model
anova(model2)                  # display seq. SS for all terms in regression

anova(model2,model1)           # perform F-test of submodel relative to model1

plot(model2)                  # obtain basic diagnostic plots for fitted model
qqnorm(resid(model2))         # including normal probability plot
qqline(resid(model2))

## Obtain Predictions and Standard Errors of Mean (Fitted) Value Estimates

xnew2 <- cbind(x1,x2,x6)
predict(model2, data.frame(x=xnew2),se.fit=T)      # obtain predictions & SEs
pointwise(predict(model2,se.fit=T))      # obtain upper & lower limits of CIs

## Perform Stepwise Regression Procedure

model0 <- lm(y ~ 1)                               # set up intercept only model
step(model0, ~ x1 + x2 + x3 + x4 + x5 + x6)  # perform stepwise regression

## Basic Procedures for One Factor (one-way) ANOVA

dissolvedat <- matrix(scan("dissolve.dat"),ncol=2,byrow=T)
diss <- dissolvedat[,2]
xtemp <- dissolvedat[,1]
temp <- factor(rep(1:6,c(3,3,3,3,3,3)))  # create factor `temp' used in aov

dissolve.tf <- data.frame(temp,diss)     # create data frame used in aov
dissolve.tf

plot.factor(dissolve.tf,main="Box Plots for Each Treatment Group")

cat("\n","One-Way ANOVA Results for Y (Diss) on factor Temp",fill=T)

aov.dissolve <- aov( diss ~ temp, dissolve.tf)   # Perform ANOVA
summary(aov.dissolve)                            # Display basic ANOVA result

cat("\n","Table of Treatment Group Means in One-Way ANOVA for Y (Diss) on factor temp",fill=T)
model.tables(aov.dissolve, type="means", se=T)

plot(fitted(aov.dissolve),resid(aov.dissolve),xlab = "Fitted Values", ylab = "Residuals", ylim = range(resid(aov.dissolve)), main="Residuals vs. Fitted Values")
abline(h=0)

qqnorm(resid(aov.dissolve),main="Normal Probability Plot of Residuals")
qqline(resid(aov.dissolve))

## Basic Procedures for Two Factor (two-way) ANOVA

q()                                               ## Ends Splus session