data <- read.table("olerZ2.csv", header = T, sep = ",")

library(MASS)

model.lower.nocov <- lm(ZAnx ~ Z84, data = data)
summary(model.lower.nocov)
model.upper.nocov <- lm(ZAnx ~ Z1 + Z2 + Z3 + Z4 + Z5 + Z6 + Z7 + Z8 + Z9 + Z10 + Z11 + Z12 + Z13 + Z14 + Z15 + Z16 + Z17 + Z18 + Z19 + Z20 + Z21 + Z22 + Z23 + Z24 + Z25 + Z26 + Z27 + Z28 + Z29 + Z30 + Z31 + Z32 + Z33 + Z34 + Z35 + Z36 + Z37 + Z38 + Z39 + Z40 + Z41 + Z42 + Z43 + Z44 + Z45 + Z46 + Z47 + Z48 + Z49 + Z50 + Z51 + Z52 + Z53 + Z54 + Z55 + Z56 + Z57 + Z58 + Z59 + Z60 + Z61 + Z62 + Z63 + Z64 + Z65 + Z66 + Z67 + Z68 + Z69 + Z70 + Z71 + Z72 + Z73 + Z74 + Z75 + Z76 + Z77 + Z78 + Z79 + Z80 + Z81 + Z82 + Z83 + Z84 + Z85 + Z86 + Z87 + Z88 + Z89 + Z90 + Z91 + Z92 + Z93 + Z94 + Z95 + Z96 + Z97 + Z98 + Z99 + Z100 + Z101 + Z102 + Z103 + Z104 + Z105 + Z106 + Z107 + Z108 + Z109 + Z110 + Z111 + Z112 + Z113 + Z114 + Z115 + Z116 + Z117 + Z118 + Z119 + Z120 + Z121 + Z122 + Z123 + Z124 + Z125 + Z126 + Z127 + Z128 + Z129 + Z130 + Z131 + Z132 + Z133 + Z134 + Z135 + Z136 + Z137 + Z138 + Z139 + Z140 + Z141 + Z142 + Z143 + Z144 + Z145 + Z146 + Z147 + Z148 + Z149 + Z150 + Z151 + Z152 + Z153 + Z154 + Z155 + Z156 + Z157, 
                      data = data)


model.lower.cov <- lm(ZAnx ~ Zsex + Zage + Z84, 
                      data = data)
model.upper.cov <- lm(ZAnx ~ Zsex + Zage + 
                             Z1 + Z2 + Z3 + Z4 + Z5 + Z6 + Z7 + Z8 + Z9 + Z10 + Z11 + Z12 + Z13 + Z14 + Z15 + Z16 + Z17 + Z18 + Z19 + Z20 + Z21 + Z22 + Z23 + Z24 + Z25 + Z26 + Z27 + Z28 + Z29 + Z30 + Z31 + Z32 + Z33 + Z34 + Z35 + Z36 + Z37 + Z38 + Z39 + Z40 + Z41 + Z42 + Z43 + Z44 + Z45 + Z46 + Z47 + Z48 + Z49 + Z50 + Z51 + Z52 + Z53 + Z54 + Z55 + Z56 + Z57 + Z58 + Z59 + Z60 + Z61 + Z62 + Z63 + Z64 + Z65 + Z66 + Z67 + Z68 + Z69 + Z70 + Z71 + Z72 + Z73 + Z74 + Z75 + Z76 + Z77 + Z78 + Z79 + Z80 + Z81 + Z82 + Z83 + Z84 + Z85 + Z86 + Z87 + Z88 + Z89 + Z90 + Z91 + Z92 + Z93 + Z94 + Z95 + Z96 + Z97 + Z98 + Z99 + Z100 + Z101 + Z102 + Z103 + Z104 + Z105 + Z106 + Z107 + Z108 + Z109 + Z110 + Z111 + Z112 + Z113 + Z114 + Z115 + Z116 + Z117 + Z118 + Z119 + Z120 + Z121 + Z122 + Z123 + Z124 + Z125 + Z126 + Z127 + Z128 + Z129 + Z130 + Z131 + Z132 + Z133 + Z134 + Z135 + Z136 + Z137 + Z138 + Z139 + Z140 + Z141 + Z142 + Z143 + Z144 + Z145 + Z146 + Z147 + Z148 + Z149 + Z150 + Z151 + Z152 + Z153 + Z154 + Z155 + Z156 + Z157, 
                      data = data)

model.upper.nocov <- lm(ZAnx ~ Z84, data = data)

step.both.nocov <- stepAIC(object = model.lower.nocov, direction = "both", scope = list(lower = model.lower.nocov, upper = model.upper.nocov))
summary(step.both.nocov)

step.both.cov <- stepAIC(object = model.lower.cov, direction = "both", scope = list(lower = model.lower.cov, upper = model.upper.cov))
summary(step.both.nocov)

step.both.nocov.klogn <- stepAIC(object = model.lower.nocov, k = log(nrow(data)), direction = "both", scope = list(lower = model.lower.nocov, upper = model.upper.nocov))
summary(step.both.nocov.klogn)

step.forward.nocov.klogn <- stepAIC(object = model.lower.nocov, k = log(nrow(data)), direction = "forward", scope = list(lower = model.lower.nocov, upper = model.upper.nocov))
summary(step.forward.nocov.klogn)

step.backward.nocov.klogn <- stepAIC(object = model.upper.nocov, k = log(nrow(data)), direction = "backward", scope = list(lower = model.lower.nocov, upper = model.upper.nocov))
summary(step.backward.nocov.klogn)


library(bootStepAIC)
bootstep.forward.nocov.klogn <- boot.stepAIC(object = model.lower.nocov, k = log(nrow(data)), direction = "forward", 
                                              scope = list(lower = model.lower.nocov, upper = model.upper.nocov),
                                              data = data, 
                                              B = 10, alpha = 0.05)
bootstep.forward.nocov.klogn$Covariates

bootstep.both.nocov.klogn <- boot.stepAIC(object = model.lower.nocov, k = log(nrow(data)), direction = "both", 
                                          scope = list(lower = model.lower.nocov, upper = model.upper.nocov),
                                          data = data, 
                                          B = 10, alpha = 0.05)
bootstep.both.nocov.klogn$Covariates


# With covariates
bootstep.forward.cov.klogn <- boot.stepAIC(object = model.lower.cov, k = log(nrow(data)), direction = "forward", 
                                           scope = list(lower = model.lower.cov, upper = model.upper.cov),
                                           data = data, 
                                           B = 1000, alpha = 0.05)
bootstep.forward.cov.klogn$Covariates

bootstep.both.cov.klogn <- boot.stepAIC(object = model.lower.cov, k = log(nrow(data)), direction = "both", 
                                        scope = list(lower = model.lower.cov, upper = model.upper.cov),
                                        data = data, 
                                        B = 1000, alpha = 0.05)
bootstep.both.cov.klogn$Covariates