davidaarmstrong

* Choosing the Optimal Confidence Level for Visual Testing
* David A. Armstrong II and William Poirier 
* 
* Requires Stata >= 16 because frames are used to save results.
* Requires the function matsort be installed with:
*     net install matsort.pkg
capture program drop testCI
program testCI
	syntax , [lev1(real .25) lev2(real .99) incr(real .01) a(real .05) easythresh(real .05) inc0 remc usemargins]
	if `lev1' <= 0  | `lev1' >= 1 {

davidaarmstrong

# Goal: Make a function that takes an interaction mdoel and transforms the output to look like subset models.
#' @parameter obj linear model object
#' @parameter data Data frame used to estimate the model
#' @parameter incl_difference Logical argument identifying whether differences between coefficient vectors should be included. 
#' @parameter ... Other arguments, currently unimplemented
interaction_to_subset <- function(obj, data, incl_difference=FALSE, ...){

}

davidaarmstrong

Y <- tibble::tribble(
  ~"V1",~"V2",~"V3",~"V4",~"V5",~"V6",~"V7",~"V8",~"V9",~"V10",~"V11",~"V12",~"V13",~"V14",~"V15",~"V16",~"V17",~"V18",~"V19",~"V20",~"V21",~"V22",~"V23",~"V24",~"V25",~"V26",~"V27", 
  209.93,209.93,209.93,209.97,209.97,211.12,211.12,211.12,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,211.15,
  201.05,201.05,201.05,201.68,201.68,204.54,204.54,204.54,205.87,209.48,209.52,209.61,209.61,209.61,210.65,210.65,210.68,210.68,210.68,210.72,210.72,210.72,210.72,210.72,210.72,210.72,210.72,
  202.72,202.72,202.72,203.2,203.2,206.38,206.38,206.38,207.6,210.31,210.35,210.44,210.44,210.44,211.5,211.5,211.54,211.54,211.54,211.58,211.58,211.58,211.58,211.58,211.58,211.58,211.58,
  202.37,202.37,202.37,202.88,202.88,206,206,206,207.24,210.13,210.17,210.27,210.27,210.27,211.33,211.33,211.36,211.36,211.4,211.4,211.4,211.4,211.4,211.4,211.4,211.4,211.4,
  202.37,202.37,202.37,202.88,202.88,206,206,206,207.24,210.13,210.17,210

davidaarmstrong

predEffects <- function(obj, digits=3){
    ## obj is a linear model object
    preds <- predict(obj, type="terms")
    rg <- apply(preds, 2, range)
    diff.rg <- apply(rg, 2, diff)
    out <- data.frame(
        min = rg[1,], 
        max = rg[2,], 
        diff = diff.rg, 
        iqr = apply(preds, 2, function(x)diff(quantile(x, c(.25,.75)))),

davidaarmstrong

library(dplyr)
library(sf)
library(rgeos)
library(rgdal)

## Download: https://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2016-eng.cfm
## choose ArcGIS (.shp) and Federal Electoral Districts (2013 Representation Order)
## read in electoral district boundaries.
eb <- read_sf("FED_CA_2_2_ENG.shp")

davidaarmstrong

ma <- function (formula, method = c("earth"), control = list(NULL), 
    ...) {
    method <- match.arg(method)
    scall <- deparse(sys.call(), width.cutoff = 200L)
    if (!is(formula, "formula")) 
        stop("formula argument in ma() needs a formula starting with ~")
    rexpr <- grepl("gamlss", sys.calls())
    for (i in length(rexpr):1) {
        position <- i
        if (rexpr[i] == TRUE)