##### General Code for calculating predictions and prediction intervals in R


### import the data


datum=read.csv(file.choose()) #imports data
head(datum) #displays the column headings of data plus first 6 rows


### Run analysis - Cannot make predictions until model has been run


results=lm(Y~X,data=datum) #runs a linear model between 'Y' and 'X' from 'datum', call results 'results
# 'Y' is the response in your excel file. BE SURE TO CHANGE THE NAME
# 'X' is the x-variable in your excel file. BE SURE TO CHANGE THE NAME
# in example, I used results=lm(Biomass~Rainfall,data=datum)
summary(results) # gives output.


##### Generating Predictions

# in this case, I'm getting predictions for a specific x-value. 
NewX=data.frame(X=x) ### where 'X' is the x-variable, and little 'x' is the value you want to predict at
###BE SURE TO CHANGE 'X' AND 'x'!!!
### For example: NewX=data.frame(Rainfall=10.0)
head(NewX) ###make sure it made the new data.frame
predictions=predict(results,NewX,interval="prediction") 
### generates your y prediction and prediction intervals
predictions ### gives the predicitons
### 'fit' is the y-value (the prediction)
### 'lwr' is the 95% lower prediction limit
### 'upr' is the 95% upper prediction limit
### prediction interval = (upr-lwr)/2


########################################################################
########## This is extra; students don't need to be able to do  #####
#########  Stuff below here #########################################
#########################################################################



##### Generate predictions for a range of x-value for purposes of plotting. 


# in this case, I'm getting predictions across the x-range
summary(datum) ###get min and max values for X
help(seq) ### help file for seq function - useful for creating sequences of values
x=seq(from=min,to=max,by=interval) 
# where 'min' is minimum value of X, 
#'max' is maximum value of x
# 'interval' is the distance between individual x values
# be sure to input values for min, max, and interval
# in example, I used x=seq(from=0,to=10,by=0.1)
NewX=data.frame(X=x) ### where 'X' is the x-variable from your excel file. BE SURE TO CHANGE THE NAME
#Example: NewX=data.frame(Rainfall=x)
head(NewX) ###make sure it made the new data.frame

predictions=predict(results,NewX,interval="prediction") 
### generates a vector of y predictions and prediction intervals

predictions ### gives the predictions

###convert 'predictions' to a data.frame
predictions=as.data.frame(predictions) 

###Adds column of X values to predictions 
predictions$X=NewX$X $ where "X" is the name of your x variable


####### plot predictions for easy viewing - NEW CODE


plot(Y~X,data=datum) 
### plot original X and Y data - replace X and Y as appriopriate

### plot predicted y value
lines(fit~X,data=predictions) 
# Where 'X' is your x variables
# 'predictions' is what you called the prediction object

### plot lower prediction interval
lines(lwr~X,data=predictions) 
# Where 'X' is your x variables
# 'predictions' is what you called the prediction object


### plot upper prediction interval
lines(upr~X,data=predictions) 
# Where 'X' is your x variables
# 'predictions' is what you called the prediction object




##### plot predictions for easy viewing - OLD CODE


plot(Y~X,data=datum,xlim=c(min,max),ylim=c(min,max)) ### plot original X and Y data - replace X and Y as appriopriate
# Be sure to specify x and y limits
# 'min' and 'max' are minimum and maximum axes values as desired
# xlab='' removes x-axis label for layered graphs
# ylab='' removes y-axis label for layered graphs
par(new=TRUE) ### ensures plot isn't replaced by next plot command; instead new info is added to plot
plot(fit~X,data=predictions,type='l',xlim=c(min,max),ylim=c(min,max),xlab='',ylab='') ### plots best fit line as a line
par(new=TRUE) ### ensures plot isn't replaced by next plot command; instead new info is added to plot
plot(lwr~X,data=predictions,type='l',xlim=c(min,max),ylim=c(min,max),xlab='',ylab='') ### plots lower prediction interval
par(new=TRUE) ### ensures plot isn't replaced by next plot command; instead new info is added to plot
plot(upr~X,data=predictions,type='l',xlim=c(min,max),ylim=c(min,max),xlab='',ylab='') ### plots upper prediction interval