For each model :
vmA1 <- sandwich::vcovHC(modelA1, type = 'HC1')

then in stargazer:

se=c(list(sqrt(diag(vmA1))),
                          list(sqrt(diag(vmA2))), ...

Suggestions

Documentation, replication etc:

• Cite properly Qc125.com as data source in README for generate_table3.R
• Move dataset variable description from the README into a codebook. Can be a separate .md file
• #4
• Write repo description in GitHub
• #3

Style:

• #5
• If there is a particular order you should run R scripts, can number them like 1_tidy_data.R

Details:

• Remove this line if not relevant anymore
  

  article_riding-volatility/codes/figure4_fragilityVSvolatility.R

  Line 41 in c70a479

  #guides(alpha = guide_legend(order = 1)) +

• Comment out this line OR move to README
  

  article_riding-volatility/codes/generate_table3.R

  Line 3 in c70a479

  devtools::install_github("clessn/clessnverse") ## if necessary

Highlights

• Looked at the JOP dataverse and this repo is quite well organized compared to some of their content!
• Very easy to read and understand what's happening without running the scripts
• Files are well organised and documented

I searched a lot and could not find a package that had functions to predict the outcome of ordinal multilevel models. So I coded the functions in mrp/functions.R with what I learned about the equations behind ordinal regressions. Would love to get some feedback on it.

I was not able to find a performant model. I am pretty sure it is due to the independent variables not being enough to predict the dependent variable, but if someone has some tips to improve them that would be great