DSJ

gistfile1.txt

#################################################################
##### Perform the paired bootstrap for linear regression    #####
#####                                                       #####
##### - yourData is a data frame containing the response and#####
#####     explanatory variables(no extra variables should be#####
#####     included)                                         #####
##### - responseName is the name of your response variable  #####
#####     as it appears in the data frame(character)        #####
##### - numberBoot is the number of bootstrap samples to    #####
#####     take                                              #####
#####                                                       #####
#####  Notes:  If you want to include categorical variables ##### 
#####             or interactions, you need to create       #####
#####             indicator variables or interaction terms  #####
#####             manually.                                 #####
#####          returns a matrix of coefficients. The order  #####
#####            of these is 1) intercept                   #####
#####            2) order of explanatory columns in data    #####
#################################################################

pairedBootstrap <- function(yourData, responseName, numberBoot){
  
  #####CREATE USEFUL VARIABLES FOR THE BOOTSTRAP METHOD
  originalSampleSize <- dim(yourData)[1]  #the data sample size
  numberExplanatory <- dim(yourData)[2]-1 #the number of explanatory variables
  
  #explanatory variable names
  explanatoryNames <- setdiff(names(yourData),responseName) 
  
  #explanatory variable formula (for the lm function)
  explanatoryFormula <- explanatoryNames[1]
  for(i in 2:length(explanatoryNames)){
    explanatoryFormula <- paste(explanatoryFormula, "+", explanatoryNames[i])
  }
  explanatoryFormula <- paste(as.character(responseName), "~", explanatoryFormula)
  
  #store coefficient estimates for each bootstrap sample in this matrix
  bootstrapCoefficients <- matrix(0, nrow=numberBoot, ncol=numberExplanatory+1)
  
  ######IMPLEMENT THE RESAMPLING PROCEDURE
  for(i in 1:numberBoot){
    #Create a vector of indices, which will make up our new sample
    samp <- sample(1:originalSampleSize, size=originalSampleSize, replace=TRUE)
    
    #save the new sample into a data frame called newSamp
    newSamp <- yourData[samp,]
    
    #run a linear regression to obtain bootstrap coef. estimates
    newReg <- lm(as.formula(explanatoryFormula), data=newSamp)
    bootstrapCoefficients[i, ] <- coef(newReg)
  }
  return(bootstrapCoefficients)
}



library(Stat2Data)
data(Perch)
attach(Perch)
pairs(~Weight+Length+Width,Perch)
reg6=lm(Weight~Length+Width,data=Perch)
par=mfrow=(c(2,2))
plot(reg6)
dataPerch=Perch
dataPerch=dataPerch[,-1]
dataPerch$LW=dataPerch$Length*dataPerch$Width
head(dataPerch)
##########
########## pairedBootstrap <- function(dataPerch, "Weight", 2000)
#########
  Boot= pairedBootstrap(dataPerch,"Weight",2000)
hist(Boot[,4],breaks=20)