Drop each Regressor step by step

Question

Is there any chance to specify the full model once and then just to drop regressors one after the other and producing a nice stargazer table with it without having to write every regression line again and again?

data <- datasets::airquality
# Treating Month and Day as crosssectional and time fixed effects

re1 <- plm(data = data, Ozone~Solar.R+Wind+Temp,
       index=c("Month", "Day"), model="within", effect="twoways") # full model 
# this is the only regression line I actually want to write.  
# The other regressors should be automatically dropped one by one in subsequent regressions.

re2 <- plm(data = data, Ozone~Wind+Temp,
       index=c("Month", "Day"), model="within", effect="twoways") # first regressor dropped

re3 <- plm(data = data, Ozone~Solar.R+Temp,
       index=c("Month", "Day"), model="within", effect="twoways") # second regressor dropped

re4 <- plm(data = data, Ozone~Solar.R+Wind,
       index=c("Month", "Day"), model="within", effect="twoways") # third regressor dropped

stargazer(re1, re2, re3, re4, 
          type = "html", 
          title = "Dropped after one another",
          out="HopeThisWorks.html")

I have looked into the step() function but this isn't quite helping as I am not aiming to drop according to significance or anything.

How do you want to choose which regressor to drop if not by significance or AIC? Stepwise regression is dodgy enough as-is, you want a even less defensible version? — Gregor Thomas, Mar 27 '18 at 14:22
Or do you want every combination? If so, then [this is a duplicate](https://stackoverflow.com/q/28606549/903061). — Gregor Thomas, Mar 27 '18 at 14:23
I just want to drop one variable at a time to proove that significance and sign of my explanatory variable of interest do not depend on any other variable being included or excluded. In my opinion this is a question of transparency. It is not about deciding which variables to keep or not! — BeSeLuFri, Mar 27 '18 at 14:24
You've got a very specific and uncommon use-case. I don't think you'll find a ready-made solution. I recommend getting the combinations of variables you want with `combn`, building the formulas with `paste`, and putting the models is a list. — Gregor Thomas, Mar 27 '18 at 14:30

score 1 · Accepted Answer · answered Mar 27 '18 at 15:07

 re1 <- plm(data = data, Ozone~Solar.R+Wind+Temp,
   index=c("Month", "Day"), model="within", effect="twoways") # full model 


 A=lapply(paste0(".~.-",c("Solar.R","Wind","Temp")),function(x)update(re1,as.formula(x)))
[[1]]

Model Formula: Ozone ~ Wind + Temp

Coefficients:
   Wind    Temp 
-2.6933  2.3735 


[[2]]

Model Formula: Ozone ~ Solar.R + Temp

Coefficients:
 Solar.R     Temp 
0.040986 2.782978 


[[3]]

Model Formula: Ozone ~ Solar.R + Wind

Coefficients:
  Solar.R      Wind 
 0.096607 -4.841992

Now to be able to access this in the global environment: use

 list2env(setNames(A,paste0("re",seq_along(A)+1)),.GlobalEnv)

 stargazer(re1, re2, re3, re4, 
      type = "html", 
      title = "Dropped after one another",
      out="HopeThisWorks.html")

score 1 · Answer 2 · answered Mar 27 '18 at 22:29

This is a more flexible alternative:

bene <- function(data = datasets::airquality,
                 depvar = "Ozone",
                 covariates = c("Wind","Temp","Solar.R")){

  funny <- function(id){
    covs <- ifelse(is.na(id),paste0(covariates,collapse = " + "),paste0(covariates[-id],collapse = " + "))

    model <- eval(parse(text=paste0("plm(data = data,",depvar," ~ ",covs,",index=c('Month', 'Day'), model='within', effect='twoways')")))

    return(model)
  }

  xx <- capture.output(stargazer::stargazer(purrr::map(c(NA,1:length(covariates)),funny),
                                            type = "html",
                                            out = paste0("results/model.html"),
                                            star.cutoffs = c(0.05,0.01,0.001),
                                            title = paste0(depvar)))

  models <- purrr::map(c(NA,1:length(covariates)),funny)
  map(models,function(x)print(summary(x)))
}

bene(data = datasets::airquality,
     depvar = "Ozone",
     covariates = c("Wind","Temp","Solar.R"))

score 0 · Answer 3 · answered Mar 27 '18 at 14:02

You could use (an adaption of) this (so far only stepwise-working) function, which shows you the stargazer output and also saves the regression table as a html file.

stepwise_model <- function(data = "datasets::airquality",
                       depvar = "Ozone",
                       covariates = c("Wind","Temp","Solar.R")){


data_df <- eval(parse(text = data)) 


models <- c()
for(q in 1:length(covariates)){

  label <- paste0("mod_",depvar,"_",q)
  models[q] <- label
  cov <- paste0(covariates[1:q],collapse = " + ")

  eval(parse(text = paste0(label," <<- plm::plm(",depvar," ~ ",cov,",data = data_df,index=c('Month', 'Day'), model='within', effect='twoways')")))

  eval(parse(text  = paste0("print(summary(",label,"))")))
}

modellist <- eval(parse(text = paste0("list(",paste0(models,collapse = ","),")")))

xx <- capture.output(stargazer::stargazer(modellist ,
                                          type = "html",
                                          out = paste0("results/paper/models/mod_",depvar,".html"),
                                          star.cutoffs = c(0.05,0.01,0.001),
                                          title = paste0(depvar)))

}

stepwise_model()

Thanks for that - But I really need code were only one variable is dropped at a time — BeSeLuFri, Mar 27 '18 at 14:16

Drop each Regressor step by step

3 Answers3