2_Social_Rank_and_Hierarchy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-eEx1m07j9pB"
   },
   "source": [
    "# CSS Lab: Social Rank and Hierarchy\n",
    "Network analysis can be used to examine social status within a community. This lab will use data from the 1995-1997 Teenage Friends and Lifestyle study [MA1997]. This study collected friendship and other data from a group of teenage students over the course of three years.\n",
    "\n",
    "## Contents\n",
    "1. [Setup](#Setup)\n",
    "2. [Social status](#Social-status)\n",
    "    1. [Minimum violation rankings](#Minimum-violation-rankings)\n",
    "3. [Directed social network](#Directed-social-network:-TFLS)\n",
    "    2. [Minimum violation rankings for TFLS](#Minimum-violation-rankings-for-TFLS)\n",
    "3. [Approximate MVRs](#Approximate-MVRs)\n",
    "    4. [Compare approximations](#Compare-approximations)\n",
    "        1. [In-Degree](#In-Degree)\n",
    "        2. [PageRank](#PageRank)\n",
    "        3. [SpringRank](#SpringRank)\n",
    "    4. [Choose an approximation](#Choose-an-approximation)\n",
    "5. [Analyze friendships](#Analyze-friendships)\n",
    "6. [Time evolution](#Time-evolution)\n",
    "    1. [Changes in social status](#Changes-in-social-status)\n",
    "    2. [Changes over multiple years](#Changes-over-multiple-years)\n",
    "    3. [Stability of social status](#Stability-of-social-status)\n",
    "    4. [Social status and communities](#Social-status-and-communities)\n",
    "7. [Behavioral data](#Behavioral-data)\n",
    "    1. [Examine individual participants](#Examine-individual-participants)\n",
    "    2. [Correlations with social status](#Correlations-with-social-status)\n",
    "    3. [Correlations with change in social status](#Correlations-with-change-in-social-status)\n",
    "    4. [Multiple regression](#Multiple-regression)\n",
    "    5. [Gender differences](#Gender-differences)\n",
    "    6. [Gender differences: multiple regression](#Gender-differences:-multiple-regression)\n",
    "    7. [Gender differences: methods comparison](#Gender-differences:-methods-comparison)\n",
    "8. [Try it yourself](#Try-it-yourself)\n",
    "9. [References](#References)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f4Iotadgj9pF"
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Qek7FkR7j9pG"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import json\n",
    "import math\n",
    "import sys\n",
    "import time\n",
    "import urllib.request\n",
    "import re\n",
    "import itertools\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "import scipy.stats as spstats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from IPython.display import display, Javascript, HTML\n",
    "import networkx as nx\n",
    "import networkx.algorithms as nxalg\n",
    "import networkx.algorithms.community as nxcom\n",
    "import networkx.readwrite as nxrw\n",
    "\n",
    "from springrank.python.SpringRank_tools import SpringRank\n",
    "from springrank.python.tools import build_graph_from_adjacency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "g54RFNuvj9pH"
   },
   "source": [
    "The following cell imports `vis.js`, the visualization library used by this lab.\n",
    "\n",
    "**Note**: If you get an error, `Javascript Error: require is not defined`, it is likely because you are running this notebook in `jupyter lab`. You need to click `Open in...` on the top right corner of the page, select `NbClassic` instead, and then rerun this notebook. This error should not happen in the classic `jupyter notebook`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NfxsMf0nj9pI"
   },
   "outputs": [],
   "source": [
    "%%javascript\n",
    "require.config({\n",
    "    paths: {\n",
    "        vis: 'https://cdnjs.cloudflare.com/ajax/libs/vis/4.21.0/vis'\n",
    "    }\n",
    "});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hzZRGkJbj9pI"
   },
   "outputs": [],
   "source": [
    "HTML('<style>{}</style>'.format(open('custom.css').read()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "id": "gpUamITEj9pI"
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def load_tfls(wave=1):\n",
    "    G = nx.DiGraph()\n",
    "    with open(\"external/TFLS/friendship-{}.csv\".format(wave)) as f:\n",
    "        # Read header\n",
    "        labels = re.split(\",\", f.readline().strip())[1:]\n",
    "        labels = [x.strip('\"') for x in labels]\n",
    "        # Read data\n",
    "        for row, row_data in enumerate(f):\n",
    "            cells = re.split(\",\", row_data.strip())\n",
    "            row_label = cells[0].strip('\"')\n",
    "            cells = cells[1:]\n",
    "            for col, cell_data in enumerate(cells):\n",
    "                if row_label == labels[col]:\n",
    "                    continue\n",
    "                if cell_data == \"1\" or cell_data == \"2\":\n",
    "                    cell_data = 1\n",
    "                    G.add_edge(row_label, labels[col], weight=cell_data, value=cell_data)\n",
    "    nx.set_node_attributes(G, dict((v, v) for v in G.nodes()), name=\"label\")\n",
    "    return G\n",
    "\n",
    "def get_colors():\n",
    "    phi = (1 + math.sqrt(5)) / 2\n",
    "    color = []\n",
    "    span = 0.5\n",
    "    low = 0.4\n",
    "    offset = 0.5\n",
    "    N = 16\n",
    "    H = 5\n",
    "    L = math.ceil(N / H)\n",
    "    for i in range(N):\n",
    "        level = (i - i % H) / H\n",
    "        t = phi * (i + 1) + offset\n",
    "        t = t - math.floor(t)\n",
    "        l = (L - level) / L\n",
    "        y = round(255 * (l * (low + span * t)))\n",
    "        b = round(255 * (l * (low + span * (1.0 - t))))\n",
    "        color.append((y, y, b))\n",
    "    return color\n",
    "\n",
    "# Standard but not colorblind-safe categorical colors\n",
    "c20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),\n",
    "    (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),\n",
    "    (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),\n",
    "    (227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),\n",
    "    (188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]\n",
    "\n",
    "def visjs_network(\n",
    "        nodes_list,\n",
    "        edges_dict,\n",
    "        node_size_field=None,\n",
    "        node_size_multiplier=1,\n",
    "        node_shadow_enabled=True,\n",
    "        node_color_border=None,\n",
    "        edge_width_field=None,\n",
    "        edge_width=None,\n",
    "        edge_arrow_to=None,\n",
    "        edge_arrow_to_scale_factor=1,\n",
    "        edge_smooth_enabled=False,\n",
    "        edge_smooth_type=\"curvedCW\",\n",
    "        graph_id=\"vis-output\"):\n",
    "\n",
    "    # Configure nodes\n",
    "    for node in nodes_list:\n",
    "        if node_size_field is not None:\n",
    "            node['size'] = node[node_size_field] * node_size_multiplier\n",
    "        node['title'] = node['label']\n",
    "        node['font'] = { 'size': 14 }\n",
    "        node['shadow'] = node_shadow_enabled\n",
    "        if node_color_border is not None:\n",
    "            node['color'] = {\n",
    "                'background': node.get('color', None),\n",
    "                'border': node_color_border\n",
    "            }\n",
    "    node_json = json.dumps(nodes_list)\n",
    "\n",
    "    # Configure edges\n",
    "    for edge in edges_dict:\n",
    "        color = edge['color']\n",
    "        edge['color'] = {\n",
    "            'color': color\n",
    "        }\n",
    "        if edge_arrow_to is not None:\n",
    "            edge['arrows'] = {\n",
    "                'to': {\n",
    "                    'enabled': edge_arrow_to,\n",
    "                    'scaleFactor': edge_arrow_to_scale_factor\n",
    "                }\n",
    "            }\n",
    "        if edge_smooth_enabled:\n",
    "            edge['smooth'] = {\n",
    "                'enabled': True,\n",
    "                'type': edge_smooth_type\n",
    "            }\n",
    "\n",
    "    edge_json = json.dumps(edges_dict)\n",
    "\n",
    "    js = Javascript(\"\"\"\n",
    "(function(element) {{\n",
    "    require(['vis'], function(vis) {{\n",
    "      // create an array with nodes\n",
    "      var nodes = new vis.DataSet({});\n",
    "      console.log(nodes)\n",
    "\n",
    "      // create an array with edges\n",
    "      var edges = new vis.DataSet({});\n",
    "      console.log(edges)\n",
    "\n",
    "      // create a network\n",
    "      var id = 'container-' + {}\n",
    "      element.append('<div style=\"width:800px;height:800px\" id=\"' + id + '\">hi</div>')\n",
    "      var container = document.getElementById(id);\n",
    "      element.append(container)\n",
    "      var data = {{\n",
    "        nodes: nodes,\n",
    "        edges: edges\n",
    "      }};\n",
    "      var options = {{\n",
    "          physics: {{\n",
    "              enabled: false\n",
    "          }},\n",
    "          interaction: {{\n",
    "              hover: true\n",
    "          }}\n",
    "      }};\n",
    "      var network = new vis.Network(container, data, options);\n",
    "    }})\n",
    "}})(element);\n",
    "    \"\"\".format(node_json, edge_json, graph_id))\n",
    "    return js\n",
    "\n",
    "def visualize_visjs(\n",
    "        G, communities=None, colors=None, safe_colors=True, default_color=\"192,192,192\",\n",
    "        node_size_field=\"node_size\", layout=\"spring\", scale=500, pos=None,\n",
    "        groups=None, smooth=False, weight=None, labels=None, title=\"\",\n",
    "        shadow=True, shape=\"dot\", node_colors=None, node_alpha=None, edge_alpha=None,\n",
    "        edge_colors=None):\n",
    "    # Defaults\n",
    "    if labels is None:\n",
    "        labels = dict()\n",
    "    if node_colors is None:\n",
    "        node_colors = dict()\n",
    "    if node_alpha is None:\n",
    "        node_alpha = dict()\n",
    "    if edge_alpha is None:\n",
    "        edge_alpha = dict()\n",
    "    if edge_colors is None:\n",
    "        edge_colors = dict()\n",
    "    # Get list of nodes and edges\n",
    "    nodes = list(G.nodes())\n",
    "    edges = list(G.edges())\n",
    "    # Change node shapes for bipartite graph\n",
    "    if groups is None:\n",
    "        node_shapes = dict()\n",
    "        node_sizes = dict()\n",
    "    else:\n",
    "        node_shapes = dict((n, \"square\") for n in groups)\n",
    "        node_sizes = dict((n, 15) for n in groups)\n",
    "        node_colors = dict((n, \"192,128,0\") for n in groups)\n",
    "    # Per-node properties\n",
    "    nodes_dict = dict((n, {\n",
    "        \"id\": n,\n",
    "        \"label\": labels.get(n, n),\n",
    "        \"node_size\": node_sizes.get(n, 5),\n",
    "        \"shape\": node_shapes.get(n, \"dot\"),\n",
    "        \"border_width\": 2,\n",
    "        \"color\": \"rgba({},{})\".format(node_colors.get(n, default_color), node_alpha.get(n, 1.0))\n",
    "        }) for n in nodes)\n",
    "    if shape == \"cat\":\n",
    "        for key, node in nodes_dict.items():\n",
    "            node[\"shape\"] = \"image\"\n",
    "            node[\"image\"] = \"cat.png\"\n",
    "    # Generate a layout for the nodes\n",
    "    edge_smooth_enabled = smooth\n",
    "    edge_width = 4\n",
    "    edge_arrow_scale = 2\n",
    "    if communities is not None and pos is None:\n",
    "        # Generate initial positions based on community\n",
    "        phi = 3.14 / len(nodes)\n",
    "        community_node = []\n",
    "        # Create list of nodes and their communities\n",
    "        for i, com in enumerate(sorted(communities, key=lambda x: len(x))):\n",
    "            for node in com:\n",
    "                community_node.append((i, node))\n",
    "        # Sort by community and\n",
    "        community_node = sorted(community_node)\n",
    "        # Generate initial position by placing communities around a circle\n",
    "        pos = dict((d[1], (math.cos(i*phi), math.sin(i*phi))) for i, d in enumerate(community_node))\n",
    "    if layout == \"circle\":\n",
    "        pos = nx.circular_layout(G, scale=scale)\n",
    "    elif layout == \"spring\":\n",
    "        try:\n",
    "            GU = G.to_undirected()\n",
    "        except AttributeError:\n",
    "            GU = G\n",
    "        pos = nx.spring_layout(GU, k=3/math.sqrt(len(nodes)), scale=scale, pos=pos)\n",
    "    # Assign position\n",
    "    for n in nodes:\n",
    "        nodes_dict[n][\"x\"] = pos[n][0]\n",
    "        nodes_dict[n][\"y\"] = pos[n][1]\n",
    "    # Calculate bounds for scaling\n",
    "    x_min = min(pos.values(), key=lambda x: x[0])[0]\n",
    "    x_max = max(pos.values(), key=lambda x: x[0])[0]\n",
    "    y_min = min(pos.values(), key=lambda x: x[1])[1]\n",
    "    y_max = max(pos.values(), key=lambda x: x[1])[1]\n",
    "    x_range = x_max - x_min\n",
    "    y_range = y_max - y_min\n",
    "    max_range = max(x_range, y_range)\n",
    "    # If we have communities, assign color based on community\n",
    "    if colors is None:\n",
    "        if safe_colors:\n",
    "            colors = [\"{},{},{}\".format(*c) for c in get_colors()]\n",
    "        else:\n",
    "            colors = [\"{},{},{}\".format(*c) for c in c20]\n",
    "    if communities is not None:\n",
    "        for i, com in enumerate(sorted(communities, key=lambda x: len(x))):\n",
    "            for node in com:\n",
    "                try:\n",
    "                    nodes_dict[node][\"color\"] = \"rgba({},{})\".format(colors[i], node_alpha.get(node, 1.0))\n",
    "                    nodes_dict[node][\"color_index\"] = i\n",
    "                except IndexError:\n",
    "                    nodes_dict[node][\"color\"] = \"rgba({},{})\".format(default_color, node_alpha.get(node, 1.0))\n",
    "    # Update color for bipartite nodes\n",
    "    for node, node_attr in nodes_dict.items():\n",
    "        if node in node_colors:\n",
    "            node_attr[\"color\"] = \"rgba({},{})\".format(node_colors.get(node, default_color), node_alpha.get(node, 1.0))\n",
    "    # Map node labels to contiguous ids\n",
    "    node_map = dict(zip(nodes,range(len(nodes))))\n",
    "    # Determine edge colors\n",
    "    edge_colors_idx = {}\n",
    "    for source, target in edges:\n",
    "        source_color = nodes_dict[source].get(\"color_index\", None)\n",
    "        target_color = nodes_dict[target].get(\"color_index\", None)\n",
    "        if source_color == target_color and source_color is not None:\n",
    "            edge_colors_idx[(source, target)] = source_color\n",
    "    for e, c in edge_colors_idx.items():\n",
    "        if c < len(colors) and e not in edge_colors:\n",
    "            edge_colors[e] = colors[c]\n",
    "    # Per-edge properties, use contiguous ids to identify nodes\n",
    "    edge_scale = math.ceil(max_range / 200)\n",
    "    edges_dict = []\n",
    "    for source, target, data in G.edges(data=True):\n",
    "        edge = {\n",
    "            \"from\": source,\n",
    "            \"to\": target,\n",
    "            \"color\": \"rgba({},{})\".format(edge_colors.get((source,target), default_color), edge_alpha.get((source,target), 0.5)),\n",
    "            \"edge_width_field\": \"value\",\n",
    "            \"value\": data.get(\"value\", 1) * edge_scale\n",
    "        }\n",
    "        edges_dict.append(edge)\n",
    "    # Convert nodes dict to node list\n",
    "    nodes_list = [nodes_dict[n] for n in nodes]\n",
    "    # Check for directed graph\n",
    "    if G.__class__ == nx.classes.digraph.DiGraph:\n",
    "        directed = True\n",
    "    else:\n",
    "        directed = False\n",
    "    # Call visjs\n",
    "    return visjs_network(\n",
    "        nodes_list, edges_dict,\n",
    "        node_size_field=\"node_size\",\n",
    "        node_size_multiplier=10.0,\n",
    "        node_shadow_enabled=shadow,\n",
    "        node_color_border=\"rgba(0,0,0,{})\".format(1.0),\n",
    "        edge_width_field=\"value\",\n",
    "        edge_width=edge_width,\n",
    "        edge_arrow_to=directed,\n",
    "        edge_arrow_to_scale_factor=edge_arrow_scale,\n",
    "        edge_smooth_enabled=edge_smooth_enabled,\n",
    "        edge_smooth_type=\"curvedCW\",\n",
    "        graph_id=hash(title))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Vg-Ijlsuj9pK"
   },
   "source": [
    "## 1. Social status\n",
    "When members of a social group can be ordered in terms of social status, it forms a pecking order. Directed network data can be used to uncover an underlying pecking order if it exists. An unreciprocated friendship can be a sign of a difference in social standing. The person who doesn't list the friendship may do so because they have too many friends to list, or because they don't want to list someone unpopular.\n",
    "\n",
    "### Minimum violation rankings\n",
    "In a perfect pecking order, friendships would only go from lower-status nodes to higher-status nodes.\n",
    "So to find a pecking order, we look for a _minimum-violation ranking_: an ordering with as few links going from high-status to low-status individuals as possible.\n",
    "Such an ordering may or may not exist.\n",
    "The more violoations in the MVR, the less heirarchical a group is.\n",
    "The example below allows you to test different orderings for the number of violations.\n",
    "\n",
    "First, visualize the network."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "x_Gc9y0bj9pL"
   },
   "source": [
    "## start from some examples with maybe only two or three people"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "id": "AhKsFX_tj9pL"
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def scale(x):\n",
    "    return (x - np.mean(x)) / np.std(x)\n",
    "\n",
    "def unit(rank):\n",
    "    v = rank.values()\n",
    "    span = max(v) - min(v)\n",
    "    return dict((k, (x - min(v)) / span) for k, x in rank.items())\n",
    "\n",
    "def giant_component(G):\n",
    "    giant_component = sorted(nxalg.weakly_connected_components(G), key=len)[0]\n",
    "    for v in set(G.nodes()) - giant_component:\n",
    "        G.remove_node(v)\n",
    "    return G\n",
    "\n",
    "def get_mvr_rank(G, time_limit=60):\n",
    "    # Iterative implementation to test all configurations\n",
    "    labels = list(G.nodes())\n",
    "    N = len(labels)\n",
    "    ranking = []\n",
    "    # Lists of remaining labels once i places have been assigned\n",
    "    remaining = [[] for i in range(N)]\n",
    "    remaining[0] = list(labels)\n",
    "    # Upper limit is all pairs\n",
    "    best_violations = N * (N - 1) / 2\n",
    "    best_ranking = []\n",
    "    start = time.time()\n",
    "    last = start\n",
    "    i = 0\n",
    "    tried = 0\n",
    "    while True:\n",
    "        i += 1\n",
    "        # Check for time limit, and sleep periodically to avoid locking cpu\n",
    "        if i % 1000 == 0:\n",
    "            t = time.time()\n",
    "            delta = t - start\n",
    "            if delta >= time_limit:\n",
    "                print(\"Time limit exceeded\")\n",
    "                print(\"{} of {:0.2e} configurations in {:0.2f}s\".format(i, math.factorial(N), delta))\n",
    "                raise RuntimeError\n",
    "            time.sleep(0)\n",
    "        try:\n",
    "            # Assign the next remaining label to the current place and move to the next\n",
    "            ranking.append(remaining[len(ranking)].pop())\n",
    "            if len(ranking) == N:\n",
    "                # Out of places, compare the full assignment to the best so far\n",
    "                trial_ranking = dict((label, i) for i, label in enumerate(ranking))\n",
    "                rank_differences, violations, reciprocated = count_violations(G, trial_ranking)\n",
    "                if violations < best_violations:\n",
    "                    best_violations = violations\n",
    "                    best_ranking = trial_ranking\n",
    "                tried += 1\n",
    "                # Backtrack\n",
    "                ranking.pop()\n",
    "            else:\n",
    "                # Move to next place\n",
    "                remaining[len(ranking)] = list(set(labels) - set(ranking))\n",
    "        except IndexError:\n",
    "            # No labels remaining, backtrack\n",
    "            try:\n",
    "                ranking.pop()\n",
    "            except IndexError:\n",
    "                # Unable to backtrack, all labels have been tried for all places\n",
    "                break\n",
    "    best_order = dict((x, i + 1) for x, i in best_ranking.items())\n",
    "    return best_ranking, best_order\n",
    "\n",
    "def get_spring_rank(G):\n",
    "    nodes = list(G.nodes())\n",
    "    A=nx.to_numpy_array(G,nodelist=nodes)\n",
    "    # Reverse spring rank so positive numbers represent higher status\n",
    "    rank = dict(zip(nodes, -1 * SpringRank(A,alpha=1.0,l0=1.0,l1=1.0)))\n",
    "    order = dict((elt[0], i) for i, elt in enumerate(sorted(rank.items(), key=lambda x: x[1], reverse = True)))\n",
    "    return rank, order\n",
    "\n",
    "def get_degree_rank(G):\n",
    "    all_nodes = set(G.nodes())\n",
    "    rank = dict((v, G.in_degree(v)) for v in all_nodes)\n",
    "    order = dict((elt[0], i) for i, elt in enumerate(sorted(rank.items(), key=lambda x: x[1], reverse = True)))\n",
    "    return rank, order\n",
    "\n",
    "def get_pagerank(G):\n",
    "    rank = nx.pagerank(G)\n",
    "    order = dict((elt[0], i) for i, elt in enumerate(sorted(rank.items(), key=lambda x: x[1], reverse = True)))\n",
    "    return rank, order\n",
    "\n",
    "def count_violations(G, order):\n",
    "    rank_differences = []\n",
    "    violations = 0\n",
    "    reciprocated = []\n",
    "    all_edges = list(G.edges())\n",
    "    for v, w in G.edges():\n",
    "        if (w, v) in all_edges:\n",
    "            reciprocated.append(True)\n",
    "        else:\n",
    "            reciprocated.append(False)\n",
    "        difference = order[v] - order[w]\n",
    "        rank_differences.append(difference)\n",
    "        if difference < 0:\n",
    "            violations += 1\n",
    "    return rank_differences, violations, reciprocated\n",
    "\n",
    "def plot_ordering(G, rank, title=\"Ordering\", scale=150, shape=\"dot\"):\n",
    "    nodes = list(G.nodes())\n",
    "    pos = dict((v, (0, r*scale)) for v, r in unit(rank).items())\n",
    "    return visualize_visjs(G, layout=None, smooth=True, pos=pos, scale=scale, title=title, shape=shape)\n",
    "\n",
    "def plot_orderings(Gs, ranks, title=\"Ordering\", connect=False, scale=150, shape=\"dot\", node_alpha=dict()):\n",
    "    pos = dict()\n",
    "    H = nx.DiGraph()\n",
    "    all_nodes = set()\n",
    "    for G in Gs:\n",
    "        all_nodes |= set(G.nodes())\n",
    "    for i, rank in enumerate(ranks):\n",
    "        rank = unit(rank)\n",
    "        GG = Gs[i]\n",
    "        GG.add_nodes_from(all_nodes)\n",
    "        new_labels = dict((v, \"{}-{}\".format(v, i)) for v in all_nodes)\n",
    "        GG = nx.relabel.relabel_nodes(GG, new_labels)\n",
    "        # Reverse y axis so higher rank is up\n",
    "        pos_i = dict((new_labels[v], (i*scale, -1*r*scale)) for v, r in rank.items())\n",
    "        pos.update(pos_i)\n",
    "        H = nx.compose(H, GG)\n",
    "    if connect:\n",
    "        for i in range(len(ranks) - 1):\n",
    "            for v in nodes:\n",
    "                H.add_edge(\"{}-{}\".format(v, i), \"{}-{}\".format(v, i+1))\n",
    "    return visualize_visjs(H, layout=None, smooth=True, pos=pos, scale=scale, title=title, shape=shape, node_alpha=node_alpha, shadow=False)\n",
    "\n",
    "def mvr_example():\n",
    "    H = nx.DiGraph()\n",
    "    H.add_edges_from([\n",
    "        ('Lemon', 'Cotton'),\n",
    "        ('Mia', 'Cotton'),\n",
    "        ('Mia', 'Lemon'),\n",
    "        ('Mia', 'Bagel'),\n",
    "        ('Lemon', 'Bagel'),\n",
    "        ('Setzer', 'Cotton'),\n",
    "        ('Lemon', 'Setzer')\n",
    "    ])\n",
    "    return H\n",
    "\n",
    "def get_communities(G):\n",
    "    communities = nx.algorithms.community.greedy_modularity_communities(G.to_undirected())\n",
    "    community = sorted(communities, key=len)\n",
    "    all_nodes = list(G.nodes())\n",
    "    df = pd.DataFrame(index=all_nodes)\n",
    "    for i, com in enumerate(communities):\n",
    "        key = \"in_community{:02d}\".format(i)\n",
    "        values = dict((v, int(v in com)) for v in all_nodes)\n",
    "        col_df = pd.DataFrame({key:values}, index=all_nodes)\n",
    "        df = pd.concat([df, col_df], axis=1)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "v8rxsh9vj9pL"
   },
   "outputs": [],
   "source": [
    "G_example = mvr_example()\n",
    "visualize_visjs(G_example, shape=\"cat\", scale = 500, title=\"MVR Example Network\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8bFlRsqvj9pM"
   },
   "source": [
    "Next, let's define an ordering, visualize it, and count the number of violations. In the visualization below, the nodes are drawn from top (highest status) to bottom (lowest status)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "92I098pDj9pM",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# To re-order the nodes, change these numbers\n",
    "mvr_order = {\n",
    "    \"Bagel\":  2,\n",
    "    \"Cotton\": 3,\n",
    "    \"Lemon\":  0,\n",
    "    \"Mia\":    1,\n",
    "    \"Setzer\": 4\n",
    "}\n",
    "# Count the violations and visualize\n",
    "rank_differences, violations, reciprocated = count_violations(G_example, mvr_order)\n",
    "print(\"Violations:\", violations)\n",
    "plot_ordering(G_example, mvr_order, scale=1000, shape=\"cat\", title=\"MVR Example\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5ju33uInj9pM"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "### Short Answer1\n",
    "In the interactive plot above, try moving the nodes around to find the fewest number of violations you can. When you've found your best order, you can change the order in the cell above and re-run the visualization.\n",
    "\n",
    "(1)In the visualization, which direction do edges point if they are violations?\n",
    "\n",
    "(2)What is the fewest number of violations you could find?\n",
    "\n",
    "(3)What is the corresponding order?\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Ib2wczdBj9pM"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KnKQ2ADYj9pN"
   },
   "source": [
    "Now let's see the actual minimum violation ranking."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7Z_FsT6Vj9pN"
   },
   "outputs": [],
   "source": [
    "mvr_score, mvr_order = get_mvr_rank(G_example)\n",
    "for label, rank in sorted(mvr_order.items(), key=lambda x: x[1]):\n",
    "    print(\"Place {}: {}\".format(rank, label))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fA76Do6ij9pN"
   },
   "source": [
    "The network above has a ranking with no violations. However, in real networks, it is uncommon to have a ranking with no violations, which is why the _minimum_ violation ranking is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MNSJHSZqj9pN"
   },
   "source": [
    "## Directed social network: TFLS\n",
    "Now we'll load and visualize friendship data from a real social network. The [MA1997] study determined friendships by asking participants to name their top friends. This method has in interesting feature: it is possible for a participant to list someone as a friend who does not list the participant as a friend. The friendship ties are _directed_. In the visualization, arrows go from the participant to the individuals they named as friends. Reciprocated friendships have arrows on both ends. The nodes and edges are also colored based on the results of a community-detection algorithm. By default, the visualization uses a colorblind-safe palette (approximately 4% of people have some form of colorblindness). If you change `safe_colors` to `False`, colors will be more distinguishable for people with full color vision.\n",
    "\n",
    "Study participants were asked to name up to 12 friends."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CtYBN1f1j9pN"
   },
   "outputs": [],
   "source": [
    "G_tfls = load_tfls()\n",
    "communities = nx.algorithms.community.greedy_modularity_communities(G_tfls.to_undirected())\n",
    "visualize_visjs(G_tfls, scale=2500, title=\"Full Network\", communities=communities, safe_colors=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WlwckM2Cj9pN"
   },
   "source": [
    "### Minimum violation ranking for TFLS\n",
    "Now let's try to find the minimum violation ranking for the Teen Friendship and Lifestyle Study. The data set is much larger than the previous example, so it may take longer to run. If it takes longer than one minute to find an MVR, the function will exit with an error message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NaxNTkUSj9pN"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    mvr_score, mvr_order = get_mvr_rank(G_tfls)\n",
    "except RuntimeError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kMRV9zx5j9pN"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "    \n",
    "### Short Answer 2\n",
    "1. Did the code above complete in less than a minute?\n",
    "2. If not, use the output of the error message to answer the following. How many possible configurations are there (note: in python, large numbers are written like 6.022e+23, which means 6.022 multiplied by 10 to the 23rd power)?\n",
    "3. How long would the algorithm have taken to complete? How does that compare to the age of the universe (4.3e+23 seconds)?\n",
    "    \n",
    "</div>   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "q5d9DgLdj9pO"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jD-dAi9yj9pO"
   },
   "source": [
    "### Approximate MVRs\n",
    "\n",
    "We'd like to find the MVR. Unfortunately, as we've seen above, finding the MVR is very difficult even with a lot of computing power (in computer science terms, it's NP-hard).\n",
    "Luckily, approximation methods exist that can run much faster. These methods find a low violation ranking but are not guaranteed to find the true minimum violation ranking.\n",
    "\n",
    "The example below explores some possible approximations. PageRank (used by Google's search engine) gives nodes higher scores for being connected to well-connected nodes. SpringRank [BLN2013] is specifically designed to find approximate MVRs. Each assigns a score to nodes, and then nodes are ranked according to that value. The visualization below shows (from left to right) degree centrality, PageRank, and SpringRank. Higher nodes have higher scores, and arrows pointing downward represent violations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ajEcxCgTj9pO"
   },
   "outputs": [],
   "source": [
    "# Calculate the rank\n",
    "dr_score, dr_order = get_degree_rank(G_tfls)\n",
    "pr_score, pr_order = get_pagerank(G_tfls)\n",
    "sr_score, sr_order = get_spring_rank(G_tfls)\n",
    "# Plot the orderings\n",
    "plot_orderings([G_tfls, G_tfls, G_tfls], [dr_score, pr_score, sr_score], title=\"MVR Examples\", scale=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Xr-61A63j9pO"
   },
   "source": [
    "### Compare approximations\n",
    "The cells below count violations for rankings based on Degree, PageRank, and SpringRank.\n",
    "\n",
    "#### There are three ways to compare these ranking methods:\n",
    "\n",
    "    number of violations: number of edges going from high-status nodes to low-status nodes.\n",
    "    \n",
    "    violation percent: the percentage of violation edges among all the edges\n",
    "    \n",
    "    mean rank difference: the average ranking differences for each pair of connected nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nim7Gdo4j9pO"
   },
   "source": [
    "#### In-Degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "uy_raNSsj9pO"
   },
   "outputs": [],
   "source": [
    "dr_differences, dr_violations, dr_reciprocated = count_violations(G_tfls, dr_order)\n",
    "print(\"Violations:\", dr_violations)\n",
    "print(\"Violation percent:\", 100 * dr_violations / len(list(G_tfls.edges())))\n",
    "print(\"Mean rank difference:\", np.mean(dr_differences))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Jpbn_kd_j9pO"
   },
   "source": [
    "#### PageRank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vhSODQGPj9pO"
   },
   "outputs": [],
   "source": [
    "pr_differences, pr_violations, pr_reciprocated = count_violations(G_tfls, pr_order)\n",
    "print(\"Violations:\", pr_violations)\n",
    "print(\"Violation percent:\", 100 * pr_violations / len(list(G_tfls.edges())))\n",
    "print(\"Mean rank difference:\", np.mean(pr_differences))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WWfgqm8lj9pO"
   },
   "source": [
    "#### SpringRank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "k7b0ubEBj9pP"
   },
   "outputs": [],
   "source": [
    "sr_differences, sr_violations, sr_reciprocated = count_violations(G_tfls, sr_order)\n",
    "print(\"Violations:\", sr_violations)\n",
    "print(\"Violation percent:\", 100 * sr_violations / len(list(G_tfls.edges())))\n",
    "print(\"Mean rank difference:\", np.mean(sr_differences))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oahVuUB6j9pP"
   },
   "source": [
    "Now let's see how highly correlated different scoring methods are. Each row/column below represents one scoring method. The off-diagonal squares plot different methods against each other, while the diagonals show histograms for each method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "e5O-DXnyj9pP"
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame({\n",
    "    \"InDegree\": dr_score,\n",
    "    \"PageRank\": pr_score,\n",
    "    \"SpringRank\": sr_score\n",
    "})\n",
    "df[\"PageRank\"] = np.log10(df[\"PageRank\"])\n",
    "pd.plotting.scatter_matrix(df, figsize=(12,12))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F3iG5_sOj9pP"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "\n",
    "### Short Answer 3\n",
    "1. Which approximation method gave the lowest number of violations?\n",
    "2. What percentage of edges violated for the best MVR approximation?\n",
    "3. If the network is very hierarchical, there should be close to 0 violations, how hierarchical is this network?\n",
    "    \n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o_UTCtnCj9pP"
   },
   "source": [
    "## 🤔 Your answer here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8KBN30Udj9pP"
   },
   "source": [
    "### Choose an approximation method\n",
    "In the cell below, select one of the scoring methods to use for the rest of the lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "p1WedFZJj9pP"
   },
   "outputs": [],
   "source": [
    "#method = \"InDegree\"\n",
    "#method = \"PageRank\"\n",
    "method = \"SpringRank\"\n",
    "\n",
    "if method == \"InDegree\":\n",
    "    order_differences, violations, reciprocated, get_rank = dr_differences, dr_violations, dr_reciprocated, get_degree_rank\n",
    "elif method == \"PageRank\":\n",
    "    order_differences, violations, reciprocated, get_rank = pr_differences, pr_violations, pr_reciprocated, get_pagerank\n",
    "elif method == \"SpringRank\":\n",
    "    order_differences, violations, reciprocated, get_rank = sr_differences, sr_violations, sr_reciprocated, get_spring_rank"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0lcIkhVTj9pT"
   },
   "source": [
    "### Analyze friendships\n",
    "\n",
    "We might wonder whether people include high or low status individuals when they list their friends. The plot below shows a histogram of the social status order differences for each friendship. Positive numbers represent someone listing a friend who has higher social status than they do. *As a reminder, high social status indicates a lower value (i.e. rank of 0 is the highest status individual).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "id": "_CpXSOXwj9pT"
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def plot_reciprocated(df):\n",
    "    rows = df[\"reciprocated\"]\n",
    "    low = min(df[\"difference\"])\n",
    "    high = max(df[\"difference\"])\n",
    "    span = high - low + 1\n",
    "    nbins = 10\n",
    "    per = span / nbins\n",
    "    x = []\n",
    "    y = []\n",
    "    err = []\n",
    "    for i in range(nbins):\n",
    "        l = round(low + i*per)\n",
    "        h = round(low + ((i+1)*per))\n",
    "        dfbin = df[(df[\"difference\"] >= l) & (df[\"difference\"] < h)]\n",
    "        N_recip = len(dfbin[dfbin[\"reciprocated\"] == True])\n",
    "        N_unrecip = len(dfbin[dfbin[\"reciprocated\"] == False])\n",
    "        x.append(round((l + h) / 2))\n",
    "        y.append(N_recip / (N_recip + N_unrecip))\n",
    "        err.append(2/math.sqrt(N_recip + N_unrecip))\n",
    "    x = np.array(x)\n",
    "    err = np.array(err)\n",
    "    f = plt.figure()\n",
    "    plt.plot(x, y)\n",
    "    plt.fill_between(x, y - err, y + err, color=\"lightblue\")\n",
    "    plt.ylim([0,1])\n",
    "    plt.tight_layout()\n",
    "    plt.xlabel(\"Friend's status difference\")\n",
    "    plt.ylabel(\"Fraction reciprocated\")\n",
    "    plt.xlim([-50, 100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "XrS8fM8sj9pT"
   },
   "outputs": [],
   "source": [
    "bins = np.arange(min(order_differences) - 0.5, max(order_differences) + 0.5)\n",
    "plt.hist(order_differences, bins=bins)\n",
    "plt.xlabel(\"Order difference\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LWlmwpTIj9pT"
   },
   "source": [
    "We are using unreciprocated friendships as an indicator of social status, <strong>so we hypothesize that an individual's friendships should be more likely to be reciprocted when they are with lower-status individuals </strong>. To confirm, the figure below shows the fraction of friendships that are reciprocated based on status difference. Positive status difference represent lower status individuals reporting friendships with higher status individuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lYOKWCCdj9pU"
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame({\"difference\": order_differences, \"reciprocated\": reciprocated})\n",
    "plot_reciprocated(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NfL94c3mj9pU"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "### Short Answer 4\n",
    "In the plot above, the line shows how often friendships are reciprocated as a function of status difference (the shaded area represents statistical uncertainty). How well does the plot confirm our hypothesis that higher status friendships are less likely to be reciprocated?\n",
    "    \n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0miaJYePj9pU"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "94KVgH41j9pU"
   },
   "source": [
    "### Time evolution\n",
    "The TFLS study surveyed participants three years in a row. These multiple measurements make it possible to analyze the change in the social network over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "id": "EWcCmJipj9pU"
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "def gradient(x):\n",
    "    y = math.floor(128 - x*64)\n",
    "    b = math.floor(128 + x*64)\n",
    "    return \"{},{},{}\".format(y, y, b)\n",
    "\n",
    "def get_deltas(Gs, scores):\n",
    "    scores = list(scores)\n",
    "    deltas = []\n",
    "    orders = []\n",
    "    for i, score in enumerate(scores):\n",
    "        # Scale to range [0,1]\n",
    "        scores[i] = unit(score)\n",
    "        # Calculate orders (small is low status)\n",
    "        orders.append(\n",
    "            dict((elt[0], i)\n",
    "            for i, elt in enumerate(sorted(scores[i].items(), key=lambda x: x[1]))))\n",
    "    # Add scores to DataFrame\n",
    "    for i in range(0, len(scores)):\n",
    "        label = []\n",
    "        delta = []\n",
    "        score = []\n",
    "        order = []\n",
    "        new_score = []\n",
    "        for v, node_score in scores[i].items():\n",
    "            # Calculate change and scale to (-1, 1)\n",
    "            label.append(v)\n",
    "            score.append(node_score)\n",
    "            order.append(orders[i][v])\n",
    "        df = pd.DataFrame({\n",
    "            \"label\":label,\n",
    "            \"score_{}\".format(i): score,\n",
    "            \"order_{}\".format(i): order})\n",
    "        df = df.set_index(\"label\")\n",
    "        deltas.append(df)\n",
    "    # Add score changes to DataFrame\n",
    "    for i in range(1, len(scores)):\n",
    "        label = []\n",
    "        delta = []\n",
    "        order_delta = []\n",
    "        for v, new_score in scores[i].items():\n",
    "            # Calculate change and scale to (-1, 1)\n",
    "            label.append(v)\n",
    "            old_score = scores[i-1][v]\n",
    "            delta.append((new_score - old_score) / 2.0)\n",
    "            old_order = orders[i-1][v]\n",
    "            new_order = orders[i][v]\n",
    "            order_delta.append(new_order - old_order)\n",
    "        df = pd.DataFrame({\n",
    "            \"label\":label,\n",
    "            \"score_delta_{}\".format(i): delta,\n",
    "            \"order_delta_{}\".format(i): order_delta})\n",
    "        df = df.set_index(\"label\")\n",
    "        deltas.append(df)\n",
    "    # Calculate similarity between old and new communities\n",
    "    all_nodes = set()\n",
    "    community_t = [nx.algorithms.community.greedy_modularity_communities(x.to_undirected()) for x in Gs]\n",
    "    node_community = [{}, {}, {}]\n",
    "    for t, communities in enumerate(community_t):\n",
    "        for i, com in enumerate(communities):\n",
    "            for v in com:\n",
    "                node_community[t][v] = i\n",
    "                all_nodes.add(v)\n",
    "    all_nodes = sorted(all_nodes)\n",
    "    for t in range(len(community_t) - 1):\n",
    "        similarity = []\n",
    "        for v in all_nodes:\n",
    "            current_com = set(community_t[t][node_community[t][v]])\n",
    "            next_com = set(community_t[t+1][node_community[t+1][v]])\n",
    "            s = len(current_com & next_com) / len(current_com | next_com)\n",
    "            similarity.append(s)\n",
    "        df = pd.DataFrame({\n",
    "            \"label\": all_nodes,\n",
    "            \"community_sim_{}\".format(t): similarity\n",
    "        }).set_index(\"label\")\n",
    "        deltas.append(df)\n",
    "    return pd.concat(deltas, axis=1)\n",
    "\n",
    "def plot_ordering_evolution(Gs, ranks, title=\"Ordering\", scale=150, shape=\"dot\"):\n",
    "    pos = dict()\n",
    "    H = nx.DiGraph()\n",
    "    all_nodes = set()\n",
    "    node_colors = {}\n",
    "    edge_colors = {}\n",
    "    edge_alpha = {}\n",
    "    for i, rank in enumerate(ranks):\n",
    "        ranks[i] = unit(rank)\n",
    "    for G in Gs:\n",
    "        all_nodes |= set(G.nodes())\n",
    "    for i, rank in enumerate(ranks):\n",
    "        GG = Gs[i].copy()\n",
    "        GG.remove_edges_from(list(GG.edges()))\n",
    "        GG.add_nodes_from(all_nodes)\n",
    "        new_labels = dict((v, \"{}-{}\".format(v, i)) for v in all_nodes)\n",
    "        GG = nx.relabel.relabel_nodes(GG, new_labels)\n",
    "        pos_i = dict((new_labels[v], (i*scale, r*scale)) for v, r in rank.items())\n",
    "        pos.update(pos_i)\n",
    "        H = nx.compose(H, GG)\n",
    "    for i in range(1, len(ranks)):\n",
    "        for v, new_rank in ranks[i].items():\n",
    "            # Calculate change and scale to (-1, 1)\n",
    "            old_label = \"{}-{}\".format(v, i-1)\n",
    "            label = \"{}-{}\".format(v, i)\n",
    "            delta = - (new_rank - ranks[i-1][v])\n",
    "            color = gradient(delta)\n",
    "            alpha = abs(new_rank - ranks[i-1][v])\n",
    "            node_colors[label] = color\n",
    "            if old_label in H and label in H:\n",
    "                H.add_edge(old_label, label)\n",
    "                edge_colors[(old_label, label)] = color\n",
    "                edge_alpha[(old_label, label)] = alpha\n",
    "    node_alpha = dict((v, 0.5) for v in H.nodes())\n",
    "    return visualize_visjs(\n",
    "        H, scale=scale, layout=None, pos=pos,\n",
    "        shape=shape, node_colors=node_colors, node_alpha=node_alpha, shadow=False,\n",
    "        edge_colors=edge_colors, edge_alpha=edge_alpha,\n",
    "        title=title)\n",
    "\n",
    "def print_evolution(label, G, score_t, order_t):\n",
    "    print(\"Participant {}\".format(label))\n",
    "    for i, H in enumerate(G):\n",
    "        orders_in = sorted([str(order_t[i][v]) for v in H.predecessors(label)])\n",
    "        orders_out = sorted([str(order_t[i][v]) for v in H.successors(label)])\n",
    "        print(\"Year {}\".format(i+1))\n",
    "        label_t = \"{}-{}\".format(label, i)\n",
    "        print(\"\\tIn-Degree: \", H.in_degree(label), \" Friends' status: {}\".format(\", \".join(orders_in)))\n",
    "        print(\"\\tOut-Degree: \", H.out_degree(label), \" Friends' status: {}\".format(\", \".join(orders_out)))\n",
    "        print(\"\\tScore: {:0.2f}\".format(score_t[i][label]))\n",
    "        print(\"\\tOrder: \", order_t[i][label])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QxhJwmpkj9pU"
   },
   "source": [
    "#### Change in social status\n",
    "The plot below shows the social status of each participant in each year, with time increasing to the right\n",
    "and status increasing **downward**.\n",
    "Arrows are drawn from each node to its place in the next year.\n",
    "Participants who move **down** in status are represented by yellow nodes/arrows, and participants who move **up** are represented by blue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m2nYYkBxj9pU"
   },
   "source": [
    "## reverse the order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jd-cM5sfj9pU"
   },
   "outputs": [],
   "source": [
    "score_t = []\n",
    "order_t = []\n",
    "G = []\n",
    "all_nodes = set()\n",
    "for wave in range(3):\n",
    "    G.append(load_tfls(wave+1))\n",
    "    all_nodes |= set(G[wave].nodes())\n",
    "for wave in range(3):\n",
    "    G[wave].add_nodes_from(all_nodes)\n",
    "    score, order = get_rank(G[wave])\n",
    "    score_t.append(score)\n",
    "    order_t.append(order)\n",
    "# Plot the orderings\n",
    "plot_ordering_evolution(G, order_t, title=\"Rank Over Time\", scale=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Y8Z2RGa9j9pV"
   },
   "source": [
    "You can use the cell below to examine how an individual's network changes over time. You can change the label parameter to examine different participants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "41vMLmMqj9pV"
   },
   "outputs": [],
   "source": [
    "print_evolution('s123', G, score_t, order_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vZXQAj37j9pV"
   },
   "source": [
    "#### Changes over multiple years\n",
    "We might ask whether individuals who move to higher status continue to do so, or do they move back towards their original position. The cell below compares the change in social status between years 1 and 2 to the change between years 2 and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qhISzDvqj9pV"
   },
   "outputs": [],
   "source": [
    "deltas = get_deltas(G, order_t)\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.plot(deltas[\"order_delta_1\"], deltas[\"order_delta_2\"], '.', markersize=20, color=\"#4F4FAF4F\")\n",
    "plt.title(\"Change in Social Status\")\n",
    "plt.xlabel(\"Year 1->2\"); plt.ylabel(\"Year 2->3\")\n",
    "plt.xticks([-150, 0, 150])\n",
    "plt.yticks([-150, 0, 150])\n",
    "plt.grid()\n",
    "# Same axis, explain status direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BOEaPEEMj9pV"
   },
   "source": [
    "#### Stability of social status\n",
    "Is social status more stable at different levels? The cell below visualizes the change in status between years 1 and 2 as a function of the year 1 status."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "H2dsi9uqj9pV"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,4))\n",
    "plt.plot(deltas[\"order_0\"], deltas[\"order_delta_1\"], '.', markersize=20, color=\"#4F4FAF4F\")\n",
    "plt.title(\"Change in Social Status\")\n",
    "plt.xlabel(\"Year 1\"); plt.ylabel(\"Year 1 -> Year 2\")\n",
    "plt.xticks([0, 160])\n",
    "plt.ylim([-160, 160])\n",
    "plt.yticks([-160, 0, 160])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Aq8Tgrdbj9pV"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "### Short Answer 5\n",
    "1. The points in the above diagram are limited to a diagonal band. If someone goes up by 160 places in year one, they can't gain any more status in year two, explain why this happens?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CcCTfK2cj9pW"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "woKTJLnVj9pW"
   },
   "source": [
    "### Behavioral data\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "id": "W8Msm8_ej9pW"
   },
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "\n",
    "def load_tfls_behavior():\n",
    "    # Load alcohol data\n",
    "    alc = pd.read_csv(\"external/TFLS/alcohol.csv\").set_index(\"Unnamed: 0\")\n",
    "    alc = alc.rename(columns={\"t1\": \"alcohol1\", \"t2\": \"alcohol2\", \"t3\": \"alcohol3\"})\n",
    "    df = alc\n",
    "    # Load age data\n",
    "    x = pd.read_csv(\"external/TFLS/age.csv\").set_index(\"Unnamed: 0\")\n",
    "    x = x.rename(columns={\"x\": \"age\"})\n",
    "    df = pd.concat([df, x], axis=1)\n",
    "    # Load sex data\n",
    "    x = pd.read_csv(\"external/TFLS/sex.csv\").set_index(\"Unnamed: 0\")\n",
    "    x = x.rename(columns={\"x\": \"female_raw\"})\n",
    "    rowna = x[\"female_raw\"].isna()\n",
    "    x[\"female\"] = x[\"female_raw\"] == 2\n",
    "    x.loc[rowna, \"female\"] = False\n",
    "    x = x.drop(columns=['female_raw'])\n",
    "    df = pd.concat([df, x], axis=1)\n",
    "    # Load cannabis data\n",
    "    x = pd.read_csv(\"external/TFLS/cannabis.csv\").set_index(\"Unnamed: 0\")\n",
    "    x = x.rename(columns={\"t1\": \"cannabis1\", \"t2\": \"cannabis2\", \"t3\": \"cannabis3\"})\n",
    "    df = pd.concat([df, x], axis=1)\n",
    "    # Load money data\n",
    "    x = pd.read_csv(\"external/TFLS/money.csv\").set_index(\"Unnamed: 0\")\n",
    "    x = x.rename(columns={\"t1\": \"money1\", \"t2\": \"money2\", \"t3\": \"money3\"})\n",
    "    df = pd.concat([df, x], axis=1)\n",
    "    # Load romantic data\n",
    "    x = pd.read_csv(\"external/TFLS/romantic.csv\").set_index(\"Unnamed: 0\")\n",
    "    x = x.rename(columns={\"t1\": \"romantic1_raw\", \"t2\": \"romantic2_raw\", \"t3\": \"romantic3_raw\"})\n",
    "    rowna = x[\"romantic1_raw\"].isna()\n",
    "    x[\"romantic1\"] = x[\"romantic1_raw\"] == 2\n",
    "    x.loc[rowna, \"romantic1\"] = False\n",
    "    rowna = x[\"romantic2_raw\"].isna()\n",
    "    x[\"romantic2\"] = x[\"romantic2_raw\"] == 2\n",
    "    x.loc[rowna, \"romantic2\"] = False\n",
    "    rowna = x[\"romantic1_raw\"].isna()\n",
    "    x[\"romantic3\"] = x[\"romantic3_raw\"] == 2\n",
    "    x.loc[rowna, \"romantic3\"] = False\n",
    "    x = x.drop(columns=['romantic1_raw', 'romantic2_raw', 'romantic3_raw'])\n",
    "    df = pd.concat([df, x], axis=1)\n",
    "    # Load tobacco data\n",
    "    x = pd.read_csv(\"external/TFLS/tobacco.csv\").set_index(\"Unnamed: 0\")\n",
    "    x = x.rename(columns={\"t1\": \"tobacco1\", \"t2\": \"tobacco2\", \"t3\": \"tobacco3\"})\n",
    "    df = pd.concat([df, x], axis=1)\n",
    "    # Set name of index\n",
    "    df[\"participant\"] = df.index\n",
    "    df = df.set_index(\"participant\")\n",
    "    return df\n",
    "\n",
    "def status_correlation(df, attributes=None, community=False):\n",
    "    if attributes == None:\n",
    "        attributes = [\"age\", \"female\", \"alcohol1\", \"cannabis1\", \"tobacco1\", \"romantic1\", \"money1\"]\n",
    "    if community:\n",
    "        com_attributes = sorted([x for x in df.columns if x[:12] == \"in_community\"])\n",
    "        attributes += com_attributes\n",
    "    status = []\n",
    "    status_p = []\n",
    "    for a in sorted(attributes):\n",
    "        isna = df[a].isna() | df[\"order_0\"].isna()\n",
    "        r, p = spstats.pearsonr(df[a][~isna], df[\"score_0\"][~isna])\n",
    "        status.append(r)\n",
    "        status_p.append(p)\n",
    "        a = a.ljust(16)\n",
    "        print(\"Status ~ {} {:0.2f}\\t(p={:0.4f})\".format(a, r, p))\n",
    "    correlations = pd.DataFrame({\n",
    "        \"attribute\": attributes,\n",
    "        \"status\": status,\n",
    "        \"status_p\": status_p,\n",
    "    }).set_index(\"attribute\")\n",
    "\n",
    "def change_correlation(df, attributes=None, community=False):\n",
    "    if attributes == None:\n",
    "        attributes = [\"age\", \"female\", \"alcohol1\", \"cannabis1\", \"tobacco1\", \"romantic1\", \"money1\"]\n",
    "    if community:\n",
    "        com_attributes = sorted([x for x in df.columns if x[:12] == \"in_community\"])\n",
    "        attributes += com_attributes\n",
    "    change = []\n",
    "    change_p = []\n",
    "    for a in sorted(attributes):\n",
    "        isna = df[a].isna() | df[\"order_delta_1\"].isna()\n",
    "        r, p = spstats.pearsonr(df[a][~isna], df['order_delta_1'][~isna])\n",
    "        change.append(r)\n",
    "        change_p.append(p)\n",
    "        a = a.ljust(16)\n",
    "        print(\"Status change ~ {} {:0.2f}\\t(p={:0.4f})\".format(a, r, p))\n",
    "    correlations = pd.DataFrame({\n",
    "        \"attribute\": attributes,\n",
    "        \"change\": change,\n",
    "        \"change_p\": change_p\n",
    "    }).set_index(\"attribute\")\n",
    "\n",
    "def participant_info(behavior, label):\n",
    "    print(\"Participant: {}, Age: {:0.2f}, Female: {}\\n\".format(\n",
    "        label, behavior.loc[label, \"age\"], behavior.loc[label, \"female\"] == 1))\n",
    "    print(\"Attribute\".ljust(16), \"Year 1\\tYear 2\\tYear3\")\n",
    "    print(\"{} {}\\t{}\\t{}\".format(\n",
    "         \"Status\".ljust(16),\n",
    "         int(behavior.loc[label, \"order_0\"]),\n",
    "         int(behavior.loc[label, \"order_1\"]),\n",
    "         int(behavior.loc[label, \"order_2\"])))\n",
    "    print(\"{} {}\\t{}\\t{}\".format(\n",
    "         \"Relationship\".ljust(16),\n",
    "         int(behavior.loc[label, \"romantic1\"]),\n",
    "         int(behavior.loc[label, \"romantic2\"]),\n",
    "         int(behavior.loc[label, \"romantic3\"])))\n",
    "    print(\"{} {:0.2f}\\t{:0.2f}\\t{:0.2f}\".format(\n",
    "         \"Pocket Money\".ljust(16),\n",
    "         behavior.loc[label, \"money1\"],\n",
    "         behavior.loc[label, \"money2\"],\n",
    "         behavior.loc[label, \"money3\"]))\n",
    "    for attribute in [\"alcohol\", \"cannabis\", \"tobacco\"]:\n",
    "        print(\"{} {:d}\\t{:d}\\t{:d}\".format(\n",
    "             attribute.ljust(16),\n",
    "             int(behavior.loc[label, attribute + \"1\"]),\n",
    "             int(behavior.loc[label, attribute + \"2\"]),\n",
    "             int(behavior.loc[label, attribute + \"3\"])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "G3hEvO24j9pW"
   },
   "source": [
    "#### Examine individual participants\n",
    "The cells below allow you to examine information about specific participants and their behavior over the course of the survey. Here is a summary of the variables available:\n",
    "\n",
    "*  Age\n",
    "*  Gender\n",
    "*  Romantic relationship status: 1 (in relationship), 0 (not)\n",
    "*  Pocket money (monthly spending money in GBP)\n",
    "*  Alcohol use: 1 (non), 2 (once or twice a year), 3 (once a month), 4 (once a week) and 5 (more than once a week)\n",
    "*  Tobacco use: 1 (non), 2 (occasional) and 3 (regular)\n",
    "*  Cannabis use: 1 (non), 2 (tried once), 3 (occasional) and 4 (regular)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BDxRcSeNj9pW"
   },
   "outputs": [],
   "source": [
    "behavior = pd.concat([load_tfls_behavior(), deltas, get_communities(G[0])], axis=1)\n",
    "participant_info(behavior, 's001')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xFzCWm6Vj9pW"
   },
   "source": [
    "The next cell shows the labels for the participants who change the most from the first wave to the second. Use the cell after that to examine the change in that participants network properties and behavior over the three years of the survey.\n",
    "* Why do individuals change in status? Do they gain or lose in-degree or out-degree?\n",
    "* What types of behaviors did they have before the change, and how did those behaviors change over time?\n",
    "* For the individuals below, can you hypothesize why their social status changed?\n",
    "* If you were conducting the study, what additional information would you collect to test your hypothesis?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tnA6Nf3oj9pW"
   },
   "outputs": [],
   "source": [
    "behavior.dropna(how=\"any\").sort_values(\"order_delta_1\").head().index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vA54VfMlj9pX"
   },
   "outputs": [],
   "source": [
    "label = \"s006\"\n",
    "print_evolution(label, G, score_t, order_t); print(\"\")\n",
    "participant_info(behavior, label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Nu1B4SF7j9pX"
   },
   "source": [
    "#### Correlations with social status\n",
    "We can check whether participant demographics or behavior correlates with social status (or change in social status). The cells below show Pearson correlation coefficients for status and several collected variables. Values between 0 and 1 represent a correlation between that variable and **lower status**, with 1 being the strongest. Values between -1 and 0 represent a correlation between that variable and **higher status**, with -1 being the strongest. The p value is also given. Large p values (typically above 0.05) are not considered statistically significant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "AAfbHX0cj9pX"
   },
   "outputs": [],
   "source": [
    "status_correlation(behavior, community=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yiIB6Ua6j9pX"
   },
   "source": [
    "#### Correlations with change in social status\n",
    "We can also see if the _change_ in social status between years is correlated with any of the demographic or behavioral variables. Positive values represent a correlation between that variable and **losing status**, and negative values represent a correlation between that variable and **gaining status**.\n",
    "## add a comment for which year of the variable is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wQFnZniNj9pX"
   },
   "outputs": [],
   "source": [
    "change_correlation(behavior, community=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d0nN8w_Dj9pX"
   },
   "source": [
    "#### Multiple regression\n",
    "The Pearson correlations above only consider the relationship between two variables at a time (e.g., social status and age). There may be correlations between independent variables, such as alcohol use and tobacco use, making it difficult to tell which is most correlated with the outcome variable. To control for these confounding factors, we can use a multiple regression analysis, taking all variables into account simultateously. Specifically, we will use an Ordinary Least Squares (OLS) regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "id": "0U2Aw0mFj9pY"
   },
   "outputs": [],
   "source": [
    "def plot_regression(df, title=None):\n",
    "    low = min(df[\"coeff\"] - 0.5*df[\"ci\"])\n",
    "    high = max(df[\"coeff\"] + 0.5*df[\"ci\"])\n",
    "    most = max(high, -low)\n",
    "    x = []\n",
    "    y = []\n",
    "    xerr = []\n",
    "    xlab = []\n",
    "    for i, factor in enumerate(df.index):\n",
    "        y.append(factor)\n",
    "        x.append(df[\"coeff\"].loc[factor])\n",
    "        xerr.append(0.5 * abs(df[\"ci\"].loc[factor]))\n",
    "    plt.errorbar(x, y, xerr=xerr, fmt='o', markersize=4, capsize=4)\n",
    "    plt.xlim([-most, most])\n",
    "    # Add center spine\n",
    "    ax = plt.gca()\n",
    "    p = ax.spines['left'].get_path()\n",
    "    s = matplotlib.spines.Spine.linear_spine(ax, \"right\")\n",
    "    s.register_axis(ax.yaxis)\n",
    "    s.set_position('zero')\n",
    "    s.set_color('gray')\n",
    "    ax.add_artist(s)\n",
    "    plt.yticks(ticks=range(len(y)), labels=y)\n",
    "    plt.xlabel(\"Effect size\")\n",
    "    if title is not None:\n",
    "        plt.title(title)\n",
    "\n",
    "def plot_behavior_ols(behavior, y_column=\"order_0\", drop=None, standardize=False, n_community=None):\n",
    "\n",
    "    df = behavior.dropna(how=\"any\").copy()\n",
    "    df['female'] = df['female'].astype(int)\n",
    "    df['romantic1'] = df['romantic1'].astype(int)\n",
    "    if standardize:\n",
    "        for factor in df.columns:\n",
    "            df[factor] = (df[factor] - df[factor].mean()) / df[factor].std()\n",
    "    attributes = [\"age\", \"female\", \"romantic1\", \"money1\", \"alcohol1\", \"cannabis1\", \"tobacco1\"]\n",
    "    if n_community:\n",
    "        com_attributes = sorted([x for x in df.columns if x[:12] == \"in_community\"])\n",
    "        # Remove last community factor to prevent collinearity\n",
    "        com_attributes = com_attributes[:n_community]\n",
    "        attributes += com_attributes\n",
    "\n",
    "    if drop is not None:\n",
    "        attributes = [it for it in attributes if it not in drop]\n",
    "\n",
    "    formula = y_column + \" ~ \" + \" + \".join(attributes)\n",
    "    print(formula)\n",
    "    model = smf.ols(formula=formula, data=df).fit()\n",
    "    df = pd.DataFrame(model.params, columns=[\"coeff\"])\n",
    "    df[\"pvalues\"] = model.pvalues\n",
    "    df[\"ci\"] = model.conf_int(alpha=0.05)[1]\n",
    "\n",
    "    title = \"OLS Regression: {}\".format(y_column)\n",
    "    return plot_regression(df, title)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WsMx_mqWj9pY"
   },
   "source": [
    "#### How to Read and Change this Plot\n",
    "Each dot represents how much a factor (labeled on the y axis) contributes to the overall change in status (labeled on the x axis). This is the \"effect size,\" also known as \"Cohen's d\". The horizontal blue lines represent the 95% confidence intervals for those estimated effect sizes. If a 95% confidence interval does not include zero (the verticle black line), then the effect is statistically significant at p<0.05. *If a blue dot is to the left of the black vertical line (a negative effect), this means that the variable is associated with higher status. If the horizontal blue lines do not cross the vertical black line for this variable, that means the variable has a statistically significant relationship with higher status.**\n",
    "\n",
    "The `y_column` parameter is set to `order_0`, the number of peers with lower status than a given student in the first survey. Other possible values include `order_delta_1`, the change in order from the first survey to the second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "j6r2RkMoj9pY"
   },
   "outputs": [],
   "source": [
    "plot_behavior_ols(behavior, y_column=\"order_0\", standardize=True, n_community=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CYOaXeDkj9pY"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "\n",
    "### Short Answer 6\n",
    "In the single variable correlation, alcohol use, cannabis use, and spending money were all significantly correlated with social status. Which of these remain significant in the mutliple regression above and which are no longer significant? Are there any variables that were not significant in the single variable correlations but became significant in the multiple regression?\n",
    "    \n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dpp-jcSVj9pY"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rj49e-UHj9pY"
   },
   "source": [
    "####  Gender Differences\n",
    "\n",
    "In an early analysis of the TFLS data set, Mitchel and Amos found gendered differences in the relationship between smoking and social status [MA1997].\n",
    "Specifically, when analyzing interviews with girls, they observed that high-status girls often discussed smoking, while low-status girls did not. However, when analyzing the transcripts of interviews with boys, they did not see a similar association.\n",
    "\n",
    "The cells below performs a similar analysis based on the numerical questionnairre responses. The first line restricts the data to female participants, while the second line restricts the data to male participants. *As a reminder, a negative correlation means the variable is associated with **higher status.**\n",
    "\n",
    "Do you find a similar result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NP0n2i5Ij9pZ"
   },
   "outputs": [],
   "source": [
    "status_correlation(behavior[behavior[\"female\"] == 1], attributes=[\"tobacco1\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "i47XgtB-j9pZ"
   },
   "outputs": [],
   "source": [
    "status_correlation(behavior[behavior[\"female\"] == 0], attributes=[\"tobacco1\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0sxwjR7Dj9pZ"
   },
   "source": [
    "<div class=\"alert-info\">\n",
    "\n",
    "\n",
    "### Short Answer 7\n",
    "#### Gender differences: multiple regression\n",
    "\n",
    "As before, any correlation with smoking may be confouded by other variables and we can control for confounds using a multiple regression. The following cell performs a multiple regression like before, but this time using only data from a single gender. Set the `female` variable to `1` to see the effects for girls or `0` to see the effects for boys. For each gender, consider these questions: Do you see an effect for smoking that is significantly different from 0? Do you see an effect for any other variables?\n",
    "    \n",
    "<div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VnaVzBySj9pZ"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-Nutp4qEj9pZ"
   },
   "outputs": [],
   "source": [
    "df = behavior[(behavior[\"female\"] == 1)].dropna(how=\"any\")\n",
    "plot_behavior_ols(df, y_column=\"order_0\", drop=[\"female\"], standardize=True, n_community=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XkD5MAtcj9pZ"
   },
   "source": [
    "<div class=\"alert-warning\">\n",
    "    \n",
    "## Reflection Questions\n",
    "\n",
    "1. Study participants were asked to name up to 12 friends. Friends who were not study participants were not recorded. Give two benefits and two drawbacks of this method for eliciting social networks? Can you think of another way to determine the friendship network of a group of people?\n",
    "    \n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BwyVca9Zj9pa"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3SRORNDtj9pa"
   },
   "source": [
    "<div class=\"alert-warning\">\n",
    "    \n",
    "2. Even if we find significant correlations between variables like status and money, it's not possible to say that a varaible causes high status rather than the other way around. How might a higher social status cause an increase in the amount of spending money?\n",
    "    \n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pCWKnVUfj9pa"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rUrMPnbpj9pa"
   },
   "source": [
    "<div class=\"alert-warning\">\n",
    "\n",
    "3. This lab has examined gender differences using two methods: Pearson correlation, and OLS multiple regression. The original authors used a qualitative analysis of interview transcripts.\n",
    "    \n",
    "    3.1 What were the similarities and differences between findings using different methods?\n",
    "    \n",
    "    3.2 Why might different methods show different findings? Can you reconcile these results?\n",
    "    \n",
    "    3.3 What are the benefits and drawbacks of each of these three types of analysis?\n",
    "    \n",
    "    <div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "N68g-Zy5j9pa"
   },
   "source": [
    "## 🤔 Your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "A0Duchmnj9pa"
   },
   "source": [
    "## References\n",
    "\n",
    "[MA1997] L. Michell, and A. Amos, \"Girls, pecking order and smoking.\" Social Science & Medicine 44(12), 1861-1869 (1997)\n",
    "\n",
    "[BLN2013] De Bacco, C., Larremore, D. B., & Moore, C. (2017). A physical model for efficient ranking in networks. arXiv preprint arXiv:1709.09002.\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}