s4gbars_R_logistic-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Jupyter notebook assumes that the R kernel for Jupyter (IRkernel) has been installed; see\n",
    "https://irkernel.github.io/installation/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# R code for logistic regression analysis of bar fractions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "\n",
    "This notebook is meant to be run within the full **s4g_barfractions** repository, including the associated data files.\n",
    "\n",
    "In addition, this notebook requires the following R packages:\n",
    "   * [survey](https://cran.r-project.org/package=survey)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "ERROR",
     "evalue": "Error in library(survey): there is no package called ‘survey’\n",
     "output_type": "error",
     "traceback": [
      "Error in library(survey): there is no package called ‘survey’\nTraceback:\n",
      "1. library(survey)"
     ]
    }
   ],
   "source": [
    "library(survey)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the following so that it points to the directory with the (text) data files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basedir <- \"./data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weighted logistic regression for Sample 1: log(M_star) alone"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Logistic regression for fraction of galaxies with bars as a function of stellar mass $\\log (M_{\\star} / M_{\\odot})$, using S4G galaxies in Sample 1 (spirals at $D \\leq 25$ Mpc) with stellar\n",
    "masses between $\\log M_{\\star} = 8.5$ and 11, with $V/V_{\\rm max}$ weighting to account for S4G angular diameter limit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data into table and then Survey-package design object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ff <- paste(basedir, \"barpresence_vs_logmstar_for_R_w25_m8.5-11.txt\", sep=\"\")\n",
    "logmstarBarWTable <- read.table(ff, header=TRUE)\n",
    "logmstarBarWDesign <- svydesign(ids=~0, data=logmstarBarWTable, weights=~weight)\n",
    "length(logmstarBarWTable$bar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standard linear logistic regression: bar fraction versus log of stellar mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logMstarWFit1 <- svyglm(bar ~ logmstar, design=logmstarBarWDesign, family=quasibinomial)\n",
    "summary(logMstarWFit1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quadratic linear logistic regression: bar fraction versus log of stellar mass + square of same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logMstarWFit2 <- svyglm(bar ~ logmstar + I(logmstar^2), design=logmstarBarWDesign, family=quasibinomial)\n",
    "summary(logMstarWFit2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of AIC values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "AIC(logMstarWFit1)\n",
    "AIC(logMstarWFit2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "747.73 - 762.586"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary\n",
    "\n",
    "Since the quadratic fit has $\\Delta$AIC $= -14.9$ relative to the linear fit, it is clearly preferred."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Weighted logistic regression for Sample 1: f(bar) vs log(M_star) and g-r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same as previous section, but now we do logistic regression versus both stellar mass and $g - r$ color, using a subsample\n",
    "of Sample 1 galaxies with color data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ff <- paste(basedir, \"barpresence_vs_logmstar-gmr_for_R_w25.txt\", sep=\"\")\n",
    "logmstargmrBarWTable <- read.table(ff, header=TRUE)\n",
    "gmrBarWDesign <- svydesign(ids=~0, data=logmstargmrBarWTable, weights=~weight)\n",
    "length(logmstargmrBarWTable$bar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear fit of $f_{\\rm bar}$ vs just $g - r$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmrWFit_gmr <- svyglm(bar ~ gmr, design=gmrBarWDesign, family=quasibinomial)\n",
    "summary(gmrWFit_gmr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit vs just logMstar for same sample: linear, then quadratic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# same sample, vs logmstar (linear) only\n",
    "gmrWFit_logmstar <- svyglm(bar ~ logmstar, design=gmrBarWDesign, family=quasibinomial)\n",
    "summary(gmrWFit_logmstar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# same sample, vs logmstar (quadratic) only\n",
    "gmrWFit_logmstar2 <- svyglm(bar ~ logmstar + I(logmstar^2), design=gmrBarWDesign, family=quasibinomial)\n",
    "summary(gmrWFit_logmstar2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finally, fit vs logMstar (quadratic) *and* g-r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmrWFit_gmrlogmstar2 <- svyglm(bar ~ logmstar + I(logmstar^2) + gmr, design=gmrBarWDesign, family=quasibinomial)\n",
    "summary(gmrWFit_gmrlogmstar2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of AIC values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "AIC(gmrWFit_gmr)\n",
    "AIC(gmrWFit_logmstar)\n",
    "AIC(gmrWFit_logmstar2)\n",
    "AIC(gmrWFit_gmrlogmstar2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary\n",
    "\n",
    "Best fit from AIC standpoint is quadratic logMstar (*without* $g - r$) -- note that its AIC is actually *lower*\n",
    "than the AIC for the quadratic logMstar + $g - r$ fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weighted logistic regression for Sample 1: f(bar) vs log(M_star) and log(f_gas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same as previous section, but now we do logistic regression versus both log of stellar mass and log of gas mass ratio $f{\\rm gas} = M_{\\rm HI} / M_{\\star}$, using a subsample\n",
    "of Sample 1 galaxies with H I data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basedir <- \"/Users/erwin/Documents/Working/Projects/Project_BarSizes/\"\n",
    "ff <- paste(basedir, \"barpresence_vs_logmstar-logfgas_for_R_w25.txt\", sep=\"\")\n",
    "logMstarfgasBarWTable <- read.table(ff, header=TRUE)\n",
    "logMstarfgasBarWDesign <- svydesign(ids=~0, data=logMstarfgasBarWTable, weights=~weight)\n",
    "length(logMstarfgasBarWTable$bar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit vs just log(f_gas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logMstarlogfgasWFit_fgas <- svyglm(bar ~ logfgas, design=logMstarfgasBarWDesign, family=quasibinomial)\n",
    "summary(logMstarlogfgasWFit_fgas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit vs just logMstar: linear, then quadratic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logMstarlogfgasWFit_logmstar <- svyglm(bar ~ logmstar, design=logMstarfgasBarWDesign, family=quasibinomial)\n",
    "summary(logMstarlogfgasWFit_logmstar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logMstarlogfgasWFit_logmstar2 <- svyglm(bar ~ logmstar + I(logmstar^2), design=logMstarfgasBarWDesign, family=quasibinomial)\n",
    "summary(logMstarlogfgasWFit_logmstar2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finally, fit vs logMstar (quadratic) *and* log(f_gas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logMstarlogfgasWFit_fgaslogmstar2 <- svyglm(bar ~ logmstar + I(logmstar^2) + logfgas, design=logMstarfgasBarWDesign, family=quasibinomial)\n",
    "summary(logMstarlogfgasWFit_fgaslogmstar2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of AIC values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "AIC(logMstarlogfgasWFit_fgas)\n",
    "AIC(logMstarlogfgasWFit_logmstar)\n",
    "AIC(logMstarlogfgasWFit_logmstar2)\n",
    "AIC(logMstarlogfgasWFit_fgaslogmstar2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary\n",
    "\n",
    "The quadratic fit using logMstar (without log f_gas) is clearly the best model."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.0"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "12px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": "3",
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}