0329heat.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TDtDifPhViwj",
        "outputId": "44214e74-2c72-490c-eea4-da0af3b7ac3e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Requirement already satisfied: rpy2 in /usr/local/lib/python3.10/dist-packages (3.4.2)\n",
            "Requirement already satisfied: cffi>=1.10.0 in /usr/local/lib/python3.10/dist-packages (from rpy2) (1.16.0)\n",
            "Requirement already satisfied: jinja2 in /usr/local/lib/python3.10/dist-packages (from rpy2) (3.1.3)\n",
            "Requirement already satisfied: pytz in /usr/local/lib/python3.10/dist-packages (from rpy2) (2023.4)\n",
            "Requirement already satisfied: tzlocal in /usr/local/lib/python3.10/dist-packages (from rpy2) (5.2)\n",
            "Requirement already satisfied: pycparser in /usr/local/lib/python3.10/dist-packages (from cffi>=1.10.0->rpy2) (2.21)\n",
            "Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from jinja2->rpy2) (2.1.5)\n",
            "\n",
            "R version 4.3.3 (2024-02-29) -- \"Angel Food Cake\"\n",
            "Copyright (C) 2024 The R Foundation for Statistical Computing\n",
            "Platform: x86_64-pc-linux-gnu (64-bit)\n",
            "\n",
            "R is free software and comes with ABSOLUTELY NO WARRANTY.\n",
            "You are welcome to redistribute it under certain conditions.\n",
            "Type 'license()' or 'licence()' for distribution details.\n",
            "\n",
            "  Natural language support but running in an English locale\n",
            "\n",
            "R is a collaborative project with many contributors.\n",
            "Type 'contributors()' for more information and\n",
            "'citation()' on how to cite R or R packages in publications.\n",
            "\n",
            "Type 'demo()' for some demos, 'help()' for on-line help, or\n",
            "'help.start()' for an HTML browser interface to help.\n",
            "Type 'q()' to quit R.\n",
            "\n",
            "> install.packages(c('sp', 'fmesher'), repos=c('http://cran.r-project.org'))\n",
            "Installing packages into ‘/usr/local/lib/R/site-library’\n",
            "(as ‘lib’ is unspecified)\n",
            "trying URL 'http://cran.r-project.org/src/contrib/sp_2.1-3.tar.gz'\n",
            "Content type 'application/x-gzip' length 1244605 bytes (1.2 MB)\n",
            "==================================================\n",
            "downloaded 1.2 MB\n",
            "\n",
            "trying URL 'http://cran.r-project.org/src/contrib/fmesher_0.1.5.tar.gz'\n",
            "Content type 'application/x-gzip' length 1001104 bytes (977 KB)\n",
            "==================================================\n",
            "downloaded 977 KB\n",
            "\n",
            "* installing *source* package ‘sp’ ...\n",
            "** package ‘sp’ successfully unpacked and MD5 sums checked\n",
            "** using staged installation\n",
            "** libs\n",
            "using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c Rcentroid.c -o Rcentroid.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c gcdist.c -o gcdist.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c init.c -o init.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c pip.c -o pip.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c pip2.c -o pip2.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c sp_xports.c -o sp_xports.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c surfaceArea.c -o surfaceArea.o\n",
            "gcc -I\"/usr/share/R/include\" -DNDEBUG       -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c zerodist.c -o zerodist.o\n",
            "gcc -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions -flto=auto -ffat-lto-objects -flto=auto -Wl,-z,relro -o sp.so Rcentroid.o gcdist.o init.o pip.o pip2.o sp_xports.o surfaceArea.o zerodist.o -L/usr/lib/R/lib -lR\n",
            "installing to /usr/local/lib/R/site-library/00LOCK-sp/00new/sp/libs\n",
            "** R\n",
            "** data\n",
            "** demo\n",
            "** inst\n",
            "** byte-compile and prepare package for lazy loading\n",
            "** help\n",
            "*** installing help indices\n",
            "** building package indices\n",
            "** installing vignettes\n",
            "** testing if installed package can be loaded from temporary location\n",
            "** checking absolute paths in shared objects and dynamic libraries\n",
            "** testing if installed package can be loaded from final location\n",
            "** testing if installed package keeps a record of temporary installation path\n",
            "* DONE (sp)\n",
            "* installing *source* package ‘fmesher’ ...\n",
            "** package ‘fmesher’ successfully unpacked and MD5 sums checked\n",
            "** using staged installation\n",
            "** libs\n",
            "using C++ compiler: ‘g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’\n",
            "using C++17\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c RcppExports.cpp -o RcppExports.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c Rcpp_interface.cc -o Rcpp_interface.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c basis.cc -o basis.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c fmesher_helpers.cc -o fmesher_helpers.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c ioutils.cc -o ioutils.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c locator.cc -o locator.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c mesh.cc -o mesh.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c meshc.cc -o meshc.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c predicates.cc -o predicates.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c trees.cc -o trees.o\n",
            "g++ -std=gnu++17 -I\"/usr/share/R/include\" -DNDEBUG -DFMESHER_WITH_R -I'/usr/local/lib/R/site-library/Rcpp/include'     -fpic  -g -O2 -ffile-prefix-map=/build/r-base-14Q6vq/r-base-4.3.3=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2  -c vector.cc -o vector.o\n",
            "g++ -std=gnu++17 -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions -flto=auto -ffat-lto-objects -flto=auto -Wl,-z,relro -o fmesher.so RcppExports.o Rcpp_interface.o basis.o fmesher_helpers.o ioutils.o locator.o mesh.o meshc.o predicates.o trees.o vector.o -L/usr/lib/R/lib -lR\n",
            "installing to /usr/local/lib/R/site-library/00LOCK-fmesher/00new/fmesher/libs\n",
            "** R\n",
            "** data\n",
            "*** moving datasets to lazyload DB\n",
            "** inst\n",
            "** byte-compile and prepare package for lazy loading\n",
            "** help\n",
            "*** installing help indices\n",
            "*** copying figures\n",
            "** building package indices\n",
            "** installing vignettes\n",
            "** testing if installed package can be loaded from temporary location\n",
            "** checking absolute paths in shared objects and dynamic libraries\n",
            "** testing if installed package can be loaded from final location\n",
            "** testing if installed package keeps a record of temporary installation path\n",
            "* DONE (fmesher)\n",
            "\n",
            "The downloaded source packages are in\n",
            "\t‘/tmp/RtmpEnm1Fb/downloaded_packages’\n",
            "> \n",
            "> \n",
            "\n",
            "R version 4.3.3 (2024-02-29) -- \"Angel Food Cake\"\n",
            "Copyright (C) 2024 The R Foundation for Statistical Computing\n",
            "Platform: x86_64-pc-linux-gnu (64-bit)\n",
            "\n",
            "R is free software and comes with ABSOLUTELY NO WARRANTY.\n",
            "You are welcome to redistribute it under certain conditions.\n",
            "Type 'license()' or 'licence()' for distribution details.\n",
            "\n",
            "  Natural language support but running in an English locale\n",
            "\n",
            "R is a collaborative project with many contributors.\n",
            "Type 'contributors()' for more information and\n",
            "'citation()' on how to cite R or R packages in publications.\n",
            "\n",
            "Type 'demo()' for some demos, 'help()' for on-line help, or\n",
            "'help.start()' for an HTML browser interface to help.\n",
            "Type 'q()' to quit R.\n",
            "\n",
            "> install.packages('INLA', repos='https://inla.r-inla-download.org/R/stable')\n",
            "Installing package into ‘/usr/local/lib/R/site-library’\n",
            "(as ‘lib’ is unspecified)\n",
            "trying URL 'https://inla.r-inla-download.org/R/stable/src/contrib/INLA_24.02.09.tar.gz'\n",
            "Content type 'application/x-gzip' length 81844375 bytes (78.1 MB)\n",
            "==================================================\n",
            "downloaded 78.1 MB\n",
            "\n",
            "* installing *binary* package ‘INLA’ ...\n",
            "* DONE (INLA)\n",
            "\n",
            "The downloaded source packages are in\n",
            "\t‘/tmp/Rtmpe1B4NE/downloaded_packages’\n",
            "> \n",
            "> \n"
          ]
        }
      ],
      "source": [
        "# 安装Python库rpy2\n",
        "!pip install rpy2\n",
        "\n",
        "# 安装R包INLA及其依赖包\n",
        "!R -e \"install.packages(c('sp', 'fmesher'), repos=c('http://cran.r-project.org'))\"\n",
        "!R -e \"install.packages('INLA', repos='https://inla.r-inla-download.org/R/stable')\""
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from rpy2 import robjects as ro\n",
        "from rpy2.robjects import pandas2ri\n",
        "from rpy2.robjects.conversion import localconverter\n",
        "\n",
        "# 激活Pandas DataFrame到R data.frame的自动转换\n",
        "pandas2ri.activate()\n",
        "\n",
        "# 加载INLA库\n",
        "ro.r('library(INLA)')\n",
        "\n",
        "def run_inla_model(formula, data, model_name):\n",
        "    with localconverter(ro.default_converter + pandas2ri.converter):\n",
        "        r_df = ro.conversion.py2rpy(data)\n",
        "\n",
        "    # 在R中创建滞后变量\n",
        "    if 'lag' in formula:\n",
        "        formula_parts = formula.split('~')\n",
        "        lags = []\n",
        "        for term in formula_parts[1].split('+'):\n",
        "            if 'lag' in term:\n",
        "                lag_term = term.strip().split('_')[-1]  # 例如 'lag1'\n",
        "                lag = int(''.join(filter(str.isdigit, lag_term)))  # 提取数字部分\n",
        "                lags.append(lag)\n",
        "        for lag in lags:\n",
        "            ro.r(f'r_df$HeatCount_lag{lag} <- c({\",\".join([\"NA\"]*lag)}, head(r_df$HeatCount, -{lag}))')\n",
        "\n",
        "    inla_call = f\"\"\"\n",
        "    result <- inla(formula = {formula}, data = r_df, family = 'gaussian',\n",
        "                   control.predictor = list(compute = TRUE))\n",
        "    \"\"\"\n",
        "    ro.r(inla_call)\n",
        "    ro.r(f'print(summary(result), digits = 3)')\n",
        "    mlik = ro.r('result$mlik')[0]\n",
        "    print(f\"模型: {model_name}\")\n",
        "    print(f\"边际对数似然: {mlik}\")\n",
        "    print()\n",
        "    return mlik\n",
        "\n",
        "# 加载数据\n",
        "df = pd.read_csv('/content/drive/My Drive/filtered_summer_data.csv')\n",
        "\n",
        "# 根据STATEname进行分组\n",
        "state_groups = df.groupby('STATEname')\n",
        "\n",
        "# 存储模型结果\n",
        "model_results = []\n",
        "\n",
        "for state, group_data in state_groups:\n",
        "    # 为每个州存储模型边际对数似然值\n",
        "    mliks = []\n",
        "    models = []\n",
        "\n",
        "    print(f\"处理州：{state}\")\n",
        "    # 运行不包含滞后变量的模型\n",
        "    base_formula = \"SentimentScore ~ HeatCount\"\n",
        "    print(\"运行不包含滞后变量的模型\")\n",
        "    base_mlik = run_inla_model(base_formula, group_data, \"基础模型\")\n",
        "    mliks.append(base_mlik)\n",
        "    models.append(\"基础模型\")\n",
        "\n",
        "    # 运行包含不同滞后天数的模型\n",
        "    for lag_days in range(1, 8):  # 比较1到7天的滞后\n",
        "        lagged_formula = \"SentimentScore ~ HeatCount\"\n",
        "        for i in range(1, lag_days + 1):\n",
        "            lagged_formula += f\" + HeatCount_lag{i}\"  # 每天滞后\n",
        "        model_name = f\"{state}_滞后{i}天模型\"\n",
        "        print(f\"运行包含滞后{lag_days}天的模型\")\n",
        "        mlik = run_inla_model(lagged_formula, group_data, model_name)\n",
        "        mliks.append(mlik)\n",
        "        models.append(model_name)\n",
        "\n",
        "    # 绘制边际对数似然值变化趋势\n",
        "    plt.figure(figsize=(12, 6))\n",
        "    plt.plot(models, mliks, marker='o')\n",
        "    plt.xticks(rotation=45)\n",
        "    plt.ylabel('边际对数似然值')\n",
        "    plt.title(f'{state} - 模型边际对数似然值变化趋势')\n",
        "    plt.show()\n",
        "\n",
        "    model_results.append((state, models, mliks))"
      ],
      "metadata": {
        "id": "15Wv1iDhi6A_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jN4PeIsGZsR-",
        "outputId": "9adaffce-d8fc-4a17-882c-4c284a5e552c"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "全国数据集的路径是/content/drive/My Drive/filtered_summer_data.csv"
      ],
      "metadata": {
        "id": "aNRXzb3jaGkO"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import rpy2.robjects as ro\n",
        "from rpy2.robjects import pandas2ri\n",
        "from rpy2.robjects.packages import importr\n",
        "\n",
        "# 自动转换Pandas DataFrame为R的data.frame\n",
        "pandas2ri.activate()\n",
        "\n",
        "# 加载R的INLA库\n",
        "inla = importr('INLA')\n",
        "base = importr('base')"
      ],
      "metadata": {
        "id": "ndaBP25VaAOg"
      },
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "import pandas as pd\n",
        "\n",
        "# 加载数据\n",
        "data_path = \"/content/drive/My Drive/filtered_summer_data.csv\"\n",
        "df = pd.read_csv(data_path)"
      ],
      "metadata": {
        "id": "ehbQMwL2aDxo"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "from rpy2 import robjects as ro\n",
        "from rpy2.robjects import pandas2ri\n",
        "from rpy2.rinterface_lib.embedded import RRuntimeError\n",
        "\n",
        "# 激活Pandas DataFrame到R data.frame的自动转换\n",
        "pandas2ri.activate()\n",
        "\n",
        "# 尝试加载INLA库,捕获任何错误\n",
        "try:\n",
        "    inla = ro.r('library(INLA)')\n",
        "except RRuntimeError as e:\n",
        "    print(f\"Failed to load INLA library in R: {e}\")\n",
        "\n",
        "def run_inla_model(formula, data):\n",
        "    \"\"\"\n",
        "    使用INLA运行贝叶斯层次模型。\n",
        "    :param formula: R风格的模型公式字符串,使用反引号括住列名。\n",
        "    :param data: 数据集,Pandas DataFrame格式。\n",
        "    \"\"\"\n",
        "    try:\n",
        "        # 将Pandas DataFrame转换为R的data.frame\n",
        "        r_df = pandas2ri.py2rpy(data)\n",
        "        # 设置R全局环境中的数据框\n",
        "        ro.globalenv['r_df'] = r_df\n",
        "        # 运行INLA模型\n",
        "        model_fit = ro.r(f\"\"\"\n",
        "            inla(formula = '{formula}', data = r_df, family = 'gaussian', control.predictor = list(compute = TRUE))\n",
        "        \"\"\")\n",
        "        # 打印模型摘要\n",
        "        print(ro.r('summary')(model_fit))\n",
        "    except RRuntimeError as e:\n",
        "        print(f\"Error running INLA model: {e}\")\n",
        "\n",
        "# 加载数据\n",
        "df = pd.read_csv('/content/drive/My Drive/filtered_summer_data.csv')\n",
        "\n",
        "# 将列名转换为小写,避免大小写问题\n",
        "df.columns = [col.lower() for col in df.columns]\n",
        "\n",
        "# 定义模型公式,使用反引号括住列名\n",
        "model_formulas = {\n",
        "    \"Model 1\": \"`sentimentscore` ~ `heatcount`\"\n",
        "}\n",
        "\n",
        "# 运行模型\n",
        "for model_name, formula in model_formulas.items():\n",
        "    print(f\"Running {model_name} with formula: {formula}\")\n",
        "    run_inla_model(formula, df)\n",
        "    print(\"\\n\" + \"=\"*50 + \"\\n\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bIttxGq5anmc",
        "outputId": "319cd5b3-3c0b-4691-d64e-62820fa22487"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: Loading required package: Matrix\n",
            "\n",
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: Loading required package: sp\n",
            "\n",
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: This is INLA_24.02.09 built 2024-02-09 03:35:28 UTC.\n",
            " - See www.r-inla.org/contact-us for how to get help.\n",
            " - List available models/likelihoods/etc with inla.list.models()\n",
            " - Use inla.doc(<NAME>) to access documentation\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Running Model 1 with formula: `sentimentscore` ~ `heatcount`\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/rpy2/robjects/pandas2ri.py:55: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.\n",
            "  for name, values in obj.iteritems():\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            " *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1 \n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : \n",
            "  object 'y...fake' not found\n",
            "The inla program failed and the maximum number of tries has been reached.\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Error running INLA model: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : \n",
            "  object 'y...fake' not found\n",
            "The inla program failed and the maximum number of tries has been reached.\n",
            "\n",
            "\n",
            "==================================================\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "这部分现在是有问题的，无法查看全国的数据，我们返回使用加州的数据进行实验"
      ],
      "metadata": {
        "id": "Za11V7U3g74E"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from rpy2 import robjects as ro\n",
        "from rpy2.robjects import pandas2ri\n",
        "from rpy2.robjects.conversion import localconverter\n",
        "\n",
        "# 激活Pandas DataFrame到R data.frame的自动转换\n",
        "pandas2ri.activate()\n",
        "\n",
        "# 加载INLA库\n",
        "ro.r('library(INLA)')\n",
        "\n",
        "def run_inla_model(formula, data, model_name):\n",
        "    with localconverter(ro.default_converter + pandas2ri.converter):\n",
        "        r_df = ro.conversion.py2rpy(data)\n",
        "\n",
        "    # 在R中创建滞后变量(如果需要)\n",
        "    if 'lag' in formula:\n",
        "        formula_parts = formula.split('~')\n",
        "        lags = []\n",
        "        for term in formula_parts[1].split('+'):\n",
        "            if 'lag' in term:\n",
        "                lag_term = term.strip().split('_')[-1]  # e.g., 'lag1'\n",
        "                lag = int(''.join(filter(str.isdigit, lag_term)))  # 提取数字部分\n",
        "                lags.append(lag)\n",
        "        for lag in lags:\n",
        "            ro.r(f'r_df$HeatCount_lag{lag} <- c({\",\".join([\"NA\"]*lag)}, head(r_df$HeatCount, -{lag}))')\n",
        "\n",
        "    inla_call = f\"\"\"\n",
        "    result <- inla(formula = {formula}, data = r_df, family = 'gaussian',\n",
        "                   control.predictor = list(compute = TRUE))\n",
        "    \"\"\"\n",
        "    ro.r(inla_call)\n",
        "    ro.r(f'print(summary(result), digits = 3)')\n",
        "    mlik = ro.r('result$mlik')[0]\n",
        "    print(f\"模型: {model_name}\")\n",
        "    print(f\"边际对数似然: {mlik}\")\n",
        "    print()\n",
        "    return mlik\n",
        "\n",
        "# 加载数据\n",
        "df = pd.read_csv('/content/TestData_California.csv')\n",
        "\n",
        "# 存储模型边际对数似然值\n",
        "mliks = []\n",
        "models = []\n",
        "\n",
        "# 运行不包含滞后变量的模型\n",
        "base_formula = \"SentimentScore ~ HeatCount\"\n",
        "print(\"运行不包含滞后变量的模型\")\n",
        "base_mlik = run_inla_model(base_formula, df, \"基础模型\")\n",
        "mliks.append(base_mlik)\n",
        "models.append(\"基础模型\")\n",
        "\n",
        "# 运行包含不同滞后天数的模型\n",
        "for lag_days in range(1, 8):  # 滞后1到7天\n",
        "    lagged_formula = \"SentimentScore ~ HeatCount\"\n",
        "    for i in range(1, lag_days + 1):\n",
        "        lagged_formula += f\" + HeatCount_lag{i}\"  # 每天滞后\n",
        "    model_name = f\"滞后{lag_days}天模型\"\n",
        "    print(f\"运行包含滞后{lag_days}天的模型\")\n",
        "    mlik = run_inla_model(lagged_formula, df, model_name)\n",
        "    mliks.append(mlik)\n",
        "    models.append(model_name)\n",
        "\n",
        "# 绘制边际对数似然值变化趋势\n",
        "plt.figure(figsize=(12, 6))\n",
        "plt.plot(models, mliks, marker='o')\n",
        "plt.xticks(rotation=45)\n",
        "plt.ylabel('边际对数似然值')\n",
        "plt.title('模型边际对数似然值变化趋势')\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "yUJQxZ1EhDVN",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 547
        },
        "outputId": "d52f617a-61a2-4d83-e2b1-ad255cbba8d3"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "运行不包含滞后变量的模型\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/rpy2/robjects/pandas2ri.py:55: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.\n",
            "  for name, values in obj.iteritems():\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            " *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1 \n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : \n",
            "  object 'SentimentScore' not found\n",
            "The inla program failed and the maximum number of tries has been reached.\n",
            "\n"
          ]
        },
        {
          "output_type": "error",
          "ename": "RRuntimeError",
          "evalue": "Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : \n  object 'SentimentScore' not found\nThe inla program failed and the maximum number of tries has been reached.\n",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mRRuntimeError\u001b[0m                             Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-8-92821d55d1fe>\u001b[0m in \u001b[0;36m<cell line: 51>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[0mbase_formula\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"SentimentScore ~ HeatCount\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"运行不包含滞后变量的模型\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 51\u001b[0;31m \u001b[0mbase_mlik\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrun_inla_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbase_formula\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"基础模型\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     52\u001b[0m \u001b[0mmliks\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbase_mlik\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     53\u001b[0m \u001b[0mmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"基础模型\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-8-92821d55d1fe>\u001b[0m in \u001b[0;36mrun_inla_model\u001b[0;34m(formula, data, model_name)\u001b[0m\n\u001b[1;32m     31\u001b[0m                    control.predictor = list(compute = TRUE))\n\u001b[1;32m     32\u001b[0m     \"\"\"\n\u001b[0;32m---> 33\u001b[0;31m     \u001b[0mro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minla_call\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     34\u001b[0m     \u001b[0mro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'print(summary(result), digits = 3)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m     \u001b[0mmlik\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'result$mlik'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/rpy2/robjects/__init__.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, string)\u001b[0m\n\u001b[1;32m    436\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    437\u001b[0m         \u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrinterface\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstring\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 438\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    439\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mconversion\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrpy2py\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mres\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    440\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/rpy2/robjects/functions.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    196\u001b[0m                 \u001b[0mv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    197\u001b[0m                 \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mr_k\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 198\u001b[0;31m         return (super(SignatureTranslatedFunction, self)\n\u001b[0m\u001b[1;32m    199\u001b[0m                 .__call__(*args, **kwargs))\n\u001b[1;32m    200\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/rpy2/robjects/functions.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    123\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    124\u001b[0m                 \u001b[0mnew_kwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconversion\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpy2rpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 125\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFunction\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnew_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mnew_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    126\u001b[0m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconversion\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrpy2py\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mres\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    127\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/rpy2/rinterface_lib/conversion.py\u001b[0m in \u001b[0;36m_\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m     43\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_cdata_res_to_rinterface\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunction\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 45\u001b[0;31m         \u001b[0mcdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     46\u001b[0m         \u001b[0;31m# TODO: test cdata is of the expected CType\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_cdata_to_rinterface\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/rpy2/rinterface.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    675\u001b[0m             )\n\u001b[1;32m    676\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0merror_occured\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 677\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0membedded\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_rinterface\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_geterrmessage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    678\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    679\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mRRuntimeError\u001b[0m: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : \n  object 'SentimentScore' not found\nThe inla program failed and the maximum number of tries has been reached.\n"
          ]
        }
      ]
    }
  ]
}