02 - analysis_plots_R.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "set.seed(1234)\n",
    "library(readr)\n",
    "library(tidyverse)\n",
    "library(modelr)\n",
    "library(bayestestR)\n",
    "library(caret)\n",
    "library(ez)\n",
    "library(apa)\n",
    "library(readr)\n",
    "library(RColorBrewer)\n",
    "library(scales)\n",
    "library(afex)\n",
    "library(ARTool)\n",
    "library(fBasics)\n",
    "library(car)\n",
    "library(grid)\n",
    "library(gridExtra)\n",
    "library(cowplot)\n",
    "library(ggsignif)\n",
    "library(tseries)\n",
    "library(Kendall)\n",
    "library(stats)\n",
    "library(ez)\n",
    "library(lme4)\n",
    "library(report)\n",
    "library(MASS)\n",
    "options(warn = -1)\n",
    "library(dplyr)\n",
    "library(emmeans)\n",
    "library(nlme)\n",
    "\n",
    "apatheme <- theme_bw() +\n",
    "  theme(\n",
    "    plot.title = element_text(hjust=0.5),\n",
    "    panel.grid.major = element_blank(),\n",
    "    panel.grid.minor = element_blank(),\n",
    "    panel.border = element_blank(),\n",
    "    axis.line = element_line(),\n",
    "    text = element_text( size = 30),\n",
    "    legend.justification = c(0,0),\n",
    "    strip.background = element_rect(fill=\"#cccccc\"),\n",
    "    axis.text.x=element_blank(),\n",
    "    axis.ticks.x=element_blank(),\n",
    "    axis.title.x=element_blank(),\n",
    "    axis.line.x = element_blank(),\n",
    "    legend.background=element_blank()\n",
    "      \n",
    "  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot_colors = c('#F8766D',\n",
    "                '#C49A00',\n",
    "                '#53B400',\n",
    "                '#00C094',\n",
    "                '#00B6EB',\n",
    "                '#A58AFF',\n",
    "                '#FB61D7')\n",
    "\n",
    "two_color = c('#04bcc4',\n",
    "            '#fc746c')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "tukey_detect<-function(dv,Tukey_crit=1.5){\n",
    "  IQR=IQR(dv,na.rm = TRUE)\n",
    "  Quant_25=quantile(dv,probs=0.25,na.rm = TRUE)\n",
    "  Quant_75=quantile(dv,probs=0.75,na.rm = TRUE)\n",
    "  upper=Quant_75+Tukey_crit*IQR\n",
    "  lower=Quant_25-Tukey_crit*IQR\n",
    "  outlier_Tukey=ifelse(dv>upper,1,ifelse(dv<lower,1,0))\n",
    "  as.numeric(paste(outlier_Tukey))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_noClose<-read_csv(\"./data_processed_noclose.csv\")\n",
    "colnames(df_noClose)<-make.names(colnames(df_noClose))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df<-read_csv(\"./data_processed.csv\")\n",
    "colnames(df)<-make.names(colnames(df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_app_usage<-read_csv(\"./app_usage.csv\")\n",
    "colnames(df_app_usage)<-make.names(colnames(df_app_usage))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_app_categories<-read_csv(\"./app_categories.csv\")\n",
    "colnames(df_app_categories)<-make.names(colnames(df_app_categories))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_entries_per_user_hours<-read_csv(\"./data_users_remaining_hours.csv\")\n",
    "colnames(df_entries_per_user_hours)<-make.names(colnames(df_entries_per_user_hours))\n",
    "\n",
    "df_entries_per_user_hours_result<-read_csv(\"./data_users_remaining_hours_result.csv\")\n",
    "colnames(df_entries_per_user_hours_result)<-make.names(colnames(df_entries_per_user_hours_result))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_entries_per_user_days<-read_csv(\"./data_users_remaining_days.csv\")\n",
    "colnames(df_entries_per_user_days)<-make.names(colnames(df_entries_per_user_days))\n",
    "\n",
    "df_entries_per_user_days_result<-read_csv(\"./data_users_remaining_days_result.csv\")\n",
    "colnames(df_entries_per_user_days_result)<-make.names(colnames(df_entries_per_user_days_result))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_entries_per_user_weeks<-read_csv(\"./data_users_remaining_weeks.csv\")\n",
    "colnames(df_entries_per_user_weeks)<-make.names(colnames(df_entries_per_user_weeks))\n",
    "\n",
    "df_entries_per_user_weeks_result<-read_csv(\"./data_users_remaining_weeks_result.csv\")\n",
    "colnames(df_entries_per_user_weeks_result)<-make.names(colnames(df_entries_per_user_weeks_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Descriptive Statistics about Users"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "how_many_sds = 2\n",
    "\n",
    "print(paste(\"Number of unique users: \", max(df_noClose$userIndex)+1))\n",
    "print(paste(\"Total number of interactions: \", length(df_noClose$userIndex)))\n",
    "print(paste(\"Avg. Interactions per user: \", length(df_noClose$resolution)/(max(df_noClose$userIndex)+1)))\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_weeks = max(weeksSinceStart)) %>%\n",
    "  summarise(avg_max_weeks = mean(max_weeks),\n",
    "              sd_max_weeks = sd(max_weeks),\n",
    "              max_max_weeks = max(max_weeks),\n",
    "              min_max_weeks = min(max_weeks))\n",
    "\n",
    "print(paste(\"Average total length (weeks): \", result$avg_max_weeks))\n",
    "print(paste(\"SD total length (weeks): \", result$sd_max_weeks))\n",
    "print(paste(\"Min total length (weeks): \", result$min_max_weeks))\n",
    "print(paste(\"Max total length (weeks): \", result$max_max_weeks))\n",
    "\n",
    "cutoff_max_weeks = result$avg_max_weeks + how_many_sds*result$sd_max_weeks\n",
    "mean_max_weeks = result$avg_max_weeks\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_weeks = max(weeksSinceStart)) %>%\n",
    "  arrange(desc(max_weeks)) %>% \n",
    "  head(5)\n",
    "\n",
    "# print(result)\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_weeks = max(weeksSinceStart)) %>%\n",
    "  arrange(desc(max_weeks)) %>% \n",
    "  tail(5)\n",
    "\n",
    "# print(result)\n",
    "\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_days = max(daysSinceStart)) %>%\n",
    "  summarise(avg_max_days = mean(max_days),\n",
    "              sd_max_days = sd(max_days))\n",
    "\n",
    "print(paste(\"Average total length (days): \", result$avg_max_days))\n",
    "print(paste(\"SD total length (days): \", result$sd_max_days))\n",
    "\n",
    "cutoff_max_days = result$avg_max_days + how_many_sds*result$sd_max_days\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_days = max(daysSinceStart)) %>%\n",
    "  arrange(desc(max_days)) %>% \n",
    "  head(5)\n",
    "\n",
    "# print(paste(\"Head:\", result))\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_days = max(minutesSinceStart)) %>%\n",
    "  arrange(desc(max_days)) %>% \n",
    "  tail(5)\n",
    "\n",
    "# print(paste(\"Tail:\", result))\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_hours = max(hoursSinceStart)) %>%\n",
    "  summarise(avg_max_hours = mean(max_hours),\n",
    "              sd_max_hours = sd(max_hours))\n",
    "\n",
    "print(paste(\"Average total length (hours): \", result$avg_max_hours))\n",
    "print(paste(\"SD total length (hours): \", result$sd_max_hours))\n",
    "\n",
    "cutoff_max_hours = result$avg_max_hours + how_many_sds*result$sd_max_hours\n",
    "\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_minutes = max(minutesSinceStart)) %>%\n",
    "  summarise(avg_max_minutes = mean(max_minutes),\n",
    "              sd_max_minutes = sd(max_minutes),\n",
    "               max_max_minutes = max(max_minutes),\n",
    "              min_max_minutes = min(max_minutes))\n",
    "\n",
    "print(paste(\"Average total length (minutes): \", result$avg_max_minutes))\n",
    "print(paste(\"SD total length (minutes): \", result$sd_max_minutes))\n",
    "print(paste(\"Min total length (minutes): \", result$min_max_minutes))\n",
    "print(paste(\"Max total length (minutes): \", result$max_max_minutes))\n",
    "\n",
    "cutoff_max_minutes = result$avg_max_minutes + how_many_sds*result$sd_max_minutes\n",
    "\n",
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_seconds = max(secondsSinceStart)) %>%\n",
    "  summarise(avg_max_seconds = mean(max_seconds),\n",
    "              sd_max_seconds = sd(max_seconds),\n",
    "               max_max_seconds = max(max_seconds),\n",
    "              min_max_seconds = min(max_seconds))\n",
    "\n",
    "print(paste(\"Average total length (seconds): \", result$avg_max_seconds))\n",
    "print(paste(\"SD total length (seconds): \", result$sd_max_seconds))\n",
    "print(paste(\"Min total length (seconds): \", result$min_max_seconds))\n",
    "print(paste(\"Max total length (seconds): \", result$max_max_seconds))\n",
    "\n",
    "print(paste(\"Cutoff minutes: \", cutoff_max_minutes))\n",
    "print(paste(\"Cutoff hours: \", cutoff_max_hours))\n",
    "print(paste(\"Cutoff days: \", cutoff_max_days))\n",
    "print(paste(\"Cutoff weeks: \", cutoff_max_weeks))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_total_effectiveness <- df_noClose %>%\n",
    "    group_by(continuedToApp) %>%\n",
    "    summarize(total_interactions = n(),\n",
    "                percentage = n() / nrow(df_noClose) * 100)\n",
    "\n",
    "df_total_effectiveness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df$lifeStatus <- factor(df$lifeStatus, levels=c('Break', 'Return', 'Dropout', 'Alive'))\n",
    "df$continuedToApp <- factor(df$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "max(result$max_weeks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "result <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_weeks = max(weeksSinceStart)) \n",
    "\n",
    "# result\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "# Plot the histogram\n",
    "max_time_plot <- ggplot(data = result, aes(x = max_weeks)) +\n",
    "  geom_histogram(binwidth = 1, fill = two_color[1], color = \"black\", alpha = 0.7) +\n",
    "  geom_vline(xintercept=cutoff_max_weeks, color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=mean_max_weeks, color='black', linetype = 'dashed') +\n",
    "  labs(x = \"Length of Recorded Usage (Weeks)\", y = \"Frequency\") +\n",
    "  scale_x_continuous(limits=c(NA,130), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,105),  expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "max_time_plot\n",
    "\n",
    "# ggsave('figures/max_time_weeks.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Users in the study by week"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "# Plot the histogram\n",
    "users_remaining_plot <- ggplot(df_entries_per_user_weeks, aes(x = weeksSinceStart, y = user_count)) +\n",
    "  geom_bar(stat = \"identity\", fill = two_color[1], color = 'black', alpha = 0.5) +\n",
    "  geom_vline(xintercept=cutoff_max_weeks, color='black', linetype = 'dashed') +\n",
    "  labs(x = \"Time Using one sec (Weeks)\", y = \"Number of Unique Users\") +\n",
    "  scale_x_continuous(limits=c(NA,80), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  # scale_y_continuous(limits=c(0,1200), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "users_remaining_plot\n",
    "\n",
    "# ggsave('figures/unique_users_weeks.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactions per user per week"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "# Plot the histogram\n",
    "interactions_per_user_days_plot <- ggplot(df_entries_per_user_days, aes(x = daysSinceStart, y = interactions_per_user)) +\n",
    "  geom_bar(stat = \"identity\", fill = two_color[1], color = 'black', alpha = 0.5) +\n",
    "  labs(x = \"Time Using one sec (Days)\", y = \"Interactions per User\") +\n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_days), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  # scale_y_continuous(limits=c(0,1200), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "interactions_per_user_days_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "R2nls <- function(nls.obj) {\n",
    "    if (class(nls.obj) != \"nls\") {\n",
    "        stop(\"The object must be of class `nls`.\")\n",
    "    }\n",
    "    da <- eval(nls.obj$data)\n",
    "    resp.name <- all.vars(summary(nls.obj)$formula)[1]\n",
    "    form <- paste(resp.name, \"~1\", sep = \"\")\n",
    "    m0 <- stats::lm(form, da)\n",
    "    an <- stats::anova(nls.obj, m0)\n",
    "    sqn <- stats::deviance(nls.obj)\n",
    "    sqe <- stats::deviance(m0)\n",
    "    r2 <- 1 - (sqn/sqe)\n",
    "    aov <- data.frame(fv = c(\"regression\", \"residuals\"),\n",
    "                      gl = c(-an$Df[2], an$Res.Df[1]),\n",
    "                      sq = c(-an$Sum[2], an$Res.Sum[1]))\n",
    "    aov$qm <- aov$sq/aov$gl\n",
    "    aov$F <- c(aov$qm[1]/aov$qm[2], NA)\n",
    "    aov$\"Pr(>F)\" <- c(1 - stats::pf(aov$F[1],\n",
    "                                    df1 = aov$gl[1],\n",
    "                                    df2 = aov$gl[2]),\n",
    "                      NA)\n",
    "    names(aov) <- c(\" \", \"Df\", \"Sum Sq\", \"Mean Sq\",\n",
    "                    \"F value\", \"Pr(>F)\")\n",
    "    return(list(anova = aov, R2 = r2))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_exp_fit <- df_entries_per_user_days[df_entries_per_user_days$daysSinceStart <= cutoff_max_days,]\n",
    "\n",
    "# Define the combined exponential function to fit\n",
    "combined_model <- function(x, A, B, C, D) {\n",
    "  A * exp(-B * x) + C * x + D\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_exp <- nls(interactions_per_user ~ combined_model(daysSinceStart, A, B, C, D),\n",
    "        data=df_exp_fit,\n",
    "        start = list(A = 10, B = .7, C = 0, D=10),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "summary(fit_exp)\n",
    "R2nls(fit_exp)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(interactions_per_user ~ A*daysSinceStart + B,\n",
    "        data=df_exp_fit,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "summary(fit_lin)\n",
    "R2nls(fit_lin)\n",
    "\n",
    "\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = df_exp_fit$daysSinceStart, y = df_exp_fit$interactions_per_user, ymin=df_exp_fit$interactions_per_user-df_exp_fit$std_error_interactions, ymax=df_exp_fit$interactions_per_user+df_exp_fit$std_error_interactions,\n",
    " Fitted_exp = predict(fit_exp, newdata = data.frame(daysSinceStart = df_exp_fit$daysSinceStart)),\n",
    " Fitted_lin = predict(fit_lin, newdata = data.frame(daysSinceStart = df_exp_fit$daysSinceStart)))\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[2], alpha = 0.5)+\n",
    "  geom_point(color = two_color[2], alpha=1) + # Data points\n",
    "  geom_line(aes(y = Fitted_exp), color = two_color[1], size=1) + # Fitted curve\n",
    "  annotate(\"text\", x=10, y=max(plot_data$y), label=expression(R^2 * \"=.913\"), parse=TRUE,vjust=1, hjust=0, size=7) +\n",
    "  labs(x = \"Time Using one sec (Days)\", y = \"App Open Attempts per User\") +\n",
    "  scale_x_continuous(limits=c(0,cutoff_max_days), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,30), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg\n",
    "\n",
    "# ggsave('figures/attempts_fit.pdf', width=10, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "m = anova(fit_lin, fit_exp)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "m[2,'Sum Sq']/sum(m['Res.Sum Sq'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_exp_fit <- df_entries_per_user_days_result[(df_entries_per_user_days_result$daysSinceStart <= cutoff_max_days) & (df_entries_per_user_days_result$continuedToApp == 'didNotOpenApp'),]\n",
    "\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "combined_model <- function(x, A, B, C, D) {\n",
    "  A * exp(-B * x) + C * x + D\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_exp <- nls(percentage ~ combined_model(daysSinceStart, A, B, C, D),\n",
    "        data=df_exp_fit,\n",
    "        start = list(A = 10, B = .7, C = 0, D=10),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "summary(fit_exp)\n",
    "R2nls(fit_exp)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(percentage ~ A*daysSinceStart + B,\n",
    "        data=df_exp_fit,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "summary(fit_lin)\n",
    "R2nls(fit_lin)\n",
    "\n",
    "\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = df_exp_fit$daysSinceStart, y = df_exp_fit$percentage, ymin=df_exp_fit$percentage-df_exp_fit$std_error_percentage, ymax=df_exp_fit$percentage+df_exp_fit$std_error_percentage,\n",
    " Fitted_exp = predict(fit_exp, newdata = data.frame(daysSinceStart = df_exp_fit$daysSinceStart)),\n",
    " Fitted_lin = predict(fit_lin, newdata = data.frame(daysSinceStart = df_exp_fit$daysSinceStart)))\n",
    "\n",
    "# Create a data frame for plotting\n",
    "# plot_data <- data.frame(x = df_exp_fit$daysSinceStart, y = df_exp_fit$percentage,\n",
    "# Fitted = predict(fit, newdata = data.frame(x = df_exp_fit$daysSinceStart)))\n",
    "\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[1], alpha = 0.5)+\n",
    "  geom_point(color = two_color[1], alpha=1) + # Data points\n",
    "  geom_line(aes(y = Fitted_exp), color = two_color[2], size=1) + # Fitted curve\n",
    "  annotate(\"text\", x=10, y=max(plot_data$y), label=expression(R^2 * \"=.840\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  labs(x = \"Time Using one sec (Days)\", y = \"Dismissed Attempts (%)\") +\n",
    "  scale_x_continuous(limits=c(0,cutoff_max_days), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,60), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg\n",
    "\n",
    "# ggsave('figures/percentage_fit.pdf', width=10, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "m = anova(fit_lin, fit_exp)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "m[2,'Sum Sq']/sum(m['Res.Sum Sq'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "# Plot the histogram\n",
    "interactions_per_user_weeks_plot <- ggplot(df_entries_per_user_weeks, aes(x = weeksSinceStart, y = interactions_per_user)) +\n",
    "  geom_bar(stat = \"identity\", fill = two_color[1], color = 'black', alpha = 0.5) +\n",
    "  labs(x = \"Time Using one sec (Weeks)\", y = \"Interactions per User\") +\n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_weeks), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  # scale_y_continuous(limits=c(0,1200), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "interactions_per_user_weeks_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Percentage of interactions where users continued to app or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the counts of each value of 'continuedToApp' for each 'weeksSinceStart'\n",
    "df_summary_weeks <- df_noClose %>%\n",
    "  count(weeksSinceStart, continuedToApp) %>%\n",
    "  group_by(weeksSinceStart) %>%\n",
    "  mutate(percentage = n / sum(n) * 100)\n",
    "\n",
    "df_summary_weeks$continuedToApp <- factor(df_summary_weeks$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Plot the stacked bar chart\n",
    "plt <- ggplot(df_summary_weeks, aes(x = weeksSinceStart)) +\n",
    "  geom_bar(stat = \"identity\",  aes(y = percentage, fill = continuedToApp)) +\n",
    "  theme_classic() +\n",
    "  labs(x = \"Time Using one sec (Weeks)\", y = \"Percentage\", fill=\"User Action\") +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_weeks), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/user_action_percentage_weeks.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the counts of each value of 'continuedToApp' for each 'weeksSinceStart'\n",
    "df_summary_days <- df_noClose %>%\n",
    "  count(daysSinceStart, continuedToApp) %>%\n",
    "  group_by(daysSinceStart) %>%\n",
    "  mutate(percentage = n / sum(n) * 100)\n",
    "\n",
    "df_summary_days$continuedToApp <- factor(df_summary_days$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Plot the stacked bar chart\n",
    "plt <- ggplot(df_summary_days, aes(x = daysSinceStart)) +\n",
    "  geom_bar(stat = \"identity\",  aes(y = percentage, fill = continuedToApp)) +\n",
    "  theme_classic() +\n",
    "  labs(x = \"Time Using one sec (Days)\", y = \"Percentage\", fill=\"User Action\") +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_days), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/user_action_percentage_days.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the counts of each value of 'continuedToApp' for each 'weeksSinceStart'\n",
    "df_summary_hours <- df_noClose %>%\n",
    "  count(hoursSinceStart, continuedToApp) %>%\n",
    "  group_by(hoursSinceStart) %>%\n",
    "  mutate(percentage = n / sum(n) * 100)\n",
    "\n",
    "df_summary_hours$continuedToApp <- factor(df_summary_hours$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Plot the stacked bar chart\n",
    "plt <- ggplot(df_summary_hours, aes(x = hoursSinceStart)) +\n",
    "  geom_bar(stat = \"identity\",  aes(y = percentage, fill = continuedToApp)) +\n",
    "  theme_classic() +\n",
    "  labs(x = \"Time Using one sec (Hours)\", y = \"Percentage\", fill=\"User Action\") +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_hours), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/user_action_percentage_hours.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# one sec Effectiveness and total interactions per user per week"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Create a dual-axis plot\n",
    "plot <- ggplot(df_entries_per_user_weeks, aes(x = weeksSinceStart)) +\n",
    "  geom_line(data=df_summary_weeks[df_summary_weeks$continuedToApp == 'didNotOpenApp',], aes(y = percentage*3, color = \"percentage\"), size = 1) +\n",
    "  geom_line(data=df_entries_per_user_weeks, aes(y = interactions_per_user, color = \"interactions_per_user\"), size = 1) +\n",
    "  scale_color_manual(values = c(two_color[2], two_color[1]), labels=c('Open Attempts per User', 'Dismissed Attempts (%)')) +\n",
    "  labs(x = \"Time Using one sec (Weeks)\", y = \"App Open Attemps per User\") +\n",
    "  theme_classic() +\n",
    "  scale_x_continuous(limits=c(0,cutoff_max_weeks), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,150), expand=c(0,0), breaks=pretty_breaks(), \n",
    "    sec.axis = sec_axis(~./3, name = \"Dismissed Attempts (%)\")) +\n",
    "  theme(text = element_text(size=25),  legend.position = \"bottom\", legend.title = element_blank(), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'))\n",
    "\n",
    "plot\n",
    "\n",
    "# ggsave('figures/attempts_and_effectiveness_weeks.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Create a dual-axis plot\n",
    "plot <- ggplot(df_entries_per_user_days, aes(x = daysSinceStart)) +\n",
    "  geom_line(data=df_summary_days[df_summary_days$continuedToApp == 'didNotOpenApp',], aes(y = percentage/3, color = \"percentage\"), size = 1) +\n",
    "  geom_line(data=df_entries_per_user_days, aes(y = interactions_per_user, color = \"interactions_per_user\"), size = 1) +\n",
    "  scale_color_manual(values = c(two_color[2], two_color[1]), labels=c('Open Attempts per User', 'Dismissed Attempts (%)')) +\n",
    "  labs(x = \"Time Using one sec (Days)\", y = \"App Open Attemps per User\") +\n",
    "  theme_classic() +\n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_days), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,30), expand=c(0,0), breaks=pretty_breaks(), \n",
    "    sec.axis = sec_axis(~.*2, name = \"Dismissed Attempts (%)\")) +\n",
    "  theme(text = element_text(size=25),  legend.position = \"bottom\", legend.title = element_blank(), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'))\n",
    "\n",
    "plot\n",
    "\n",
    "# ggsave('figures/attempts_and_effectiveness_days.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# df_summary_per_user_hours$percentage = df_summary_hours[df_summary_hours$continuedToApp == 'didNotOpenApp',]$percentage\n",
    "\n",
    "# df_effectiveness$continuedToApp <- factor(df_effectiveness$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "# Create a dual-axis plot\n",
    "plot <- ggplot(df_entries_per_user_hours, aes(x = hoursSinceStart)) +\n",
    "  geom_line(data=df_summary_hours[df_summary_hours$continuedToApp == 'didNotOpenApp',], aes(y = percentage/20, color = \"percentage\"), size = 1) +\n",
    "  geom_line(data=df_entries_per_user_hours, aes(y = interactions_per_user, color = \"interactions_per_user\"), size = 1) +\n",
    "  scale_color_manual(values = c(two_color[2], two_color[1]), labels=c('Open Attempts per User', 'Dismissed Attempts (%)')) +\n",
    "  labs(x = \"Time Using one sec (Hours)\", y = \"App Open Attemps per User\") +\n",
    "  theme_classic() +\n",
    "  scale_x_continuous(limits=c(NA,cutoff_max_hours), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,5), expand=c(0,0), breaks=pretty_breaks(), \n",
    "    sec.axis = sec_axis(~.*20, name = \"Dismissed Attempts (%)\")) +\n",
    "  theme(text = element_text(size=25),  legend.position = \"bottom\", legend.title = element_blank(), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'))\n",
    "\n",
    "plot\n",
    "\n",
    "# ggsave('figures/attempts_and_effectiveness_days.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_entries_per_user_weeks_result$continuedToApp <- factor(df_entries_per_user_weeks_result$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "plt <- ggplot(df_entries_per_user_weeks_result, aes(x=weeksSinceStart)) +\n",
    "geom_bar(stat = \"identity\",  aes(y = interactions_per_user, fill = continuedToApp)) +\n",
    "theme_classic() +\n",
    "labs(x = \"Time Using one sec (Weeks)\", y = \"Total App Open Attemps per User\", fill=\"User Action\") +\n",
    "scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "scale_color_manual(values = two_color) + \n",
    "scale_x_continuous(limits=c(NA,15), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "# scale_y_continuous(limits=c(0,150), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=12, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/opens_effectiveness_per_user_weeks.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_entries_per_user_days_result$continuedToApp <- factor(df_entries_per_user_days_result$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "plt <- ggplot(df_entries_per_user_days_result, aes(x=daysSinceStart)) +\n",
    "geom_bar(stat = \"identity\", aes(y = interactions_per_user, fill = continuedToApp)) +\n",
    "theme_classic() +\n",
    "labs(x = \"Time Using one sec (Days)\", y = \"Total App Open Attemps per User\", fill=\"User Action\") +\n",
    "scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "scale_color_manual(values = two_color) + \n",
    "scale_x_continuous(limits=c(NA,50), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "# scale_y_continuous(limits=c(0,30), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=12, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/opens_effectiveness_per_user_days.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_entries_per_user_hours_result$continuedToApp <- factor(df_entries_per_user_hours_result$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "plt <- ggplot(df_entries_per_user_hours_result, aes(x=hoursSinceStart)) +\n",
    "geom_bar(stat = \"identity\",  aes(y = interactions_per_user, fill = continuedToApp)) +\n",
    "theme_classic() +\n",
    "labs(x = \"Time Using one sec (Hours)\", y = \"Total App Open Attemps per User\", fill=\"User Action\") +\n",
    "scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "scale_color_manual(values = two_color) + \n",
    "scale_x_continuous(limits=c(48,96), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "# scale_y_continuous(limits=c(0,2000), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=12, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/opens_effectiveness_per_user_hours.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trends and Correlations for Attempts and Effectiveness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "trend_days <- df_entries_per_user_days[df_entries_per_user_days$daysSinceStart <= cutoff_max_days,]\n",
    "trend_percentage_days <- df_entries_per_user_days_result[(df_entries_per_user_days_result$daysSinceStart <= cutoff_max_days) & (df_entries_per_user_days_result$continuedToApp == 'didNotOpenApp'),]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(adf.test(trend_days$interactions_per_user))\n",
    "print(adf.test(trend_percentage_days$percentage))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(MannKendall(trend_days$interactions_per_user))\n",
    "print(MannKendall(trend_percentage_days$percentage))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "pearsonTest(trend_days$daysSinceStart, trend_days$interactions_per_user)\n",
    "pearsonTest(trend_percentage_days$daysSinceStart, trend_percentage_days$percentage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_break <- df_noClose[df_noClose$lifeStatus %in% c(\"Break\"), ]\n",
    "# df_dropout <- df_noClose[df_noClose$lifeStatus %in% c(\"Dropout\"), ]\n",
    "df_return <- df_noClose[df_noClose$lifeStatus %in% c(\"Return\"), ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(paste(\"Number of Breaks: \", length(df_break$userIndex)))\n",
    "\n",
    "unique_count <- df_break %>%\n",
    "  summarise(unique_count = n_distinct(userIndex))\n",
    "print(paste(\"Number of Per User (for users with at least one break): \", length(df_break$userIndex)/unique_count))\n",
    "\n",
    "\n",
    "temp <- df_break %>% \n",
    "    group_by(userIndex) %>%\n",
    "  summarise(break_count = sum(lifeStatus == \"Break\", na.rm = TRUE))\n",
    "print(paste(\"Max breaks: \", max(temp$break_count)))\n",
    "\n",
    "\n",
    "print(paste(\"Min break length (hours): \", min(df_break$hourDiffForward)))\n",
    "print(paste(\"Min break length (days): \", min(df_break$dayDiffForward)))\n",
    "print(paste(\"Mean break length (hours): \", mean(df_break$hourDiffForward)))\n",
    "print(paste(\"Mean break length (days): \", mean(df_break$dayDiffForward)))\n",
    "print(paste(\"Max break length (hours): \", max(df_break$hourDiffForward)))\n",
    "print(paste(\"Max break length (days): \", max(df_break$dayDiffForward)))\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First and Last Week"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# # Sort the dataframe by 'userIndex' and 'weeksSinceStudySetup'\n",
    "df <- df %>%\n",
    "  arrange(userIndex, realWeeksDifference)\n",
    "  \n",
    "# Get the firstWeek and lastWeek for each participant\n",
    "grouped_data <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(firstWeek = first(realWeeksDifference), lastWeek = last(realWeeksDifference))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# # Sort the dataframe by 'userIndex' and 'weeksSinceStudySetup'\n",
    "df <- df %>%\n",
    "  arrange(userIndex, realWeeksDifference)\n",
    "  \n",
    "# Get the firstWeek and lastWeek for each participant\n",
    "grouped_data <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(firstWeek = first(timestamp), lastWeek = last(timestamp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_data <- grouped_data %>%\n",
    "  arrange(firstWeek)\n",
    "  \n",
    "\n",
    "oldest_data_date <- first(grouped_data$firstWeek)\n",
    "\n",
    "print(paste('Earliest data point', oldest_data_date))\n",
    "# last(grouped_data$firstWeek)\n",
    "\n",
    "grouped_data <- grouped_data %>%\n",
    "  arrange(lastWeek)\n",
    "\n",
    "first_join_date <- first(grouped_data$lastWeek)\n",
    "last_join_date <- last(grouped_data$lastWeek)\n",
    "\n",
    "print(paste('Earliest participant join date', first_join_date))\n",
    "print(paste('Last participant join date', last_join_date))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "lims <- as.POSIXct(strptime(c(\"2020-09-01 00:00\",\"2023-05-01 00:00\"), format = \"%Y-%m-%d %H:%M\"))  \n",
    "\n",
    "# Plotting the horizontal lines with red dots at endpoints\n",
    "plt <- ggplot(grouped_data, aes(x = firstWeek, xend = lastWeek, y = userIndex, yend = userIndex)) +\n",
    "  geom_segment(size = 1, color='darkgrey') +\n",
    "  geom_point(aes(x = firstWeek, color = \"First Week\"), size = 3) +\n",
    "  geom_point(aes(x = lastWeek, color=\"Last Week\"), size = 3) +\n",
    "  # geom_vline(xintercept=as.POSIXct(\"2022-11-05 00:00:00\",tz=Sys.timezone()), color='red', linetype = 'dashed', size=1) +\n",
    "  # geom_vline(xintercept=as.POSIXct(\"2023-08-27 00:00:00\",tz=Sys.timezone()), color='black', linetype = 'dashed', size=1) +\n",
    "  geom_vline(xintercept=first_join_date, color='black', linetype = 'dashed', size=1) +\n",
    "  geom_vline(xintercept=last_join_date, color='black', linetype = 'dashed', size=1) +\n",
    "  scale_color_manual(labels=c(\"First Entry\", \"Last Entry\"), values=two_color) + \n",
    "  labs(x=\"Date\", y=\"User Index\", color=\"\")+\n",
    "  scale_x_datetime(breaks=date_breaks('3 month'), limits=lims, labels=date_format(\"%Y-%m\")) +\n",
    "  scale_y_continuous(expand=c(0,0), breaks=pretty_breaks(), limits=c(0,max(grouped_data$userIndex)+10)) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\", plot.margin = margin(t=15, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "\n",
    "\n",
    "plt\n",
    "\n",
    "# ggsave('figures/timeline_first_last.pdf', width=20, height=7)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Breaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Function to get breakStart and breakEnd values\n",
    "get_breaks <- function(weeks) {\n",
    "  diffs <- diff(weeks)\n",
    "  break_start <- weeks[-length(weeks)][diffs > 1]\n",
    "  break_end <- weeks[-1][diffs > 1]\n",
    "  return(data.frame(breakStart = break_start, breakEnd = break_end))\n",
    "}\n",
    "\n",
    "# Group by 'userIndex' and apply the get_breaks function\n",
    "breaks_df <- df_noClose %>%\n",
    "  group_by(userIndex) %>%\n",
    "  do(get_breaks(.$timestamp)) %>%\n",
    "  ungroup()\n",
    "\n",
    "# Remove rows where breakStart and breakEnd are missing (i.e., single data points for each userIndex)\n",
    "breaks_df <- breaks_df[complete.cases(breaks_df), ]\n",
    "\n",
    "# Print the resulting dataframe\n",
    "print(breaks_df)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From weeksSinceStudyStart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# # Sort the dataframe by 'userIndex' and 'weeksSinceStudySetup'\n",
    "df <- df %>%\n",
    "  arrange(userIndex, daysSinceStart)\n",
    "  \n",
    "# Get the firstWeek and lastWeek for each participant\n",
    "grouped_data <- df_noClose %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(firstDay = first(daysSinceStart), lastDay = last(daysSinceStart))\n",
    "\n",
    "\n",
    "# Sort the dataframe by 'firstWeek'\n",
    "sorted_df <- grouped_data[order(grouped_data$lastDay), ]\n",
    "\n",
    "# Create a new column 'currentIndex' containing the sorted indices\n",
    "sorted_df$currentIndex <- seq_len(nrow(sorted_df))\n",
    "\n",
    "# Plotting the horizontal lines with red dots at endpoints\n",
    "ggplot(sorted_df, aes(x = firstDay, xend = lastDay, y = currentIndex, yend = currentIndex)) +\n",
    "  geom_segment(size = 1) +\n",
    "  geom_point(aes(x = lastDay), size = 3, color=two_color[1]) +\n",
    "  # scale_color_manual(values = two_color[1]) + \n",
    "  ylab(\"userIndex\") +\n",
    "  xlab(\"Time (Days)\") +\n",
    "  ggtitle(\"Participant Timeline\") +\n",
    "  theme_classic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Function to get breakStart and breakEnd values\n",
    "get_breaks <- function(days) {\n",
    "  diffs <- diff(days)\n",
    "  break_start <- days[-length(days)][diffs > 2.7]\n",
    "  break_end <- days[-1][diffs > 2.7]\n",
    "  return(data.frame(breakStart = break_start, breakEnd = break_end))\n",
    "}\n",
    "\n",
    "# Group by 'userIndex' and apply the get_breaks function\n",
    "breaks_df <- df_noClose %>%\n",
    "  group_by(userIndex) %>%\n",
    "  do(get_breaks(.$daysSinceStart)) %>%\n",
    "  ungroup()\n",
    "\n",
    "# Remove rows where breakStart and breakEnd are missing (i.e., single data points for each userIndex)\n",
    "breaks_df <- breaks_df[complete.cases(breaks_df), ]\n",
    "\n",
    "# Print the resulting dataframe\n",
    "print(breaks_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "breaks_df$breakLength <- breaks_df$breakEnd - breaks_df$breakStart\n",
    "head(breaks_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Break Lengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# head(df_break)\n",
    "# Create the histogram using ggplot2\n",
    "\n",
    "histogram <- ggplot(data = df_break, aes(x = dayDiffForward)) +\n",
    "  geom_histogram(binwidth = 1, fill = two_color[1], color = \"black\", alpha = 0.7) +\n",
    "  labs(x = \"Break Length (Days)\", y = \"Frequency\") +\n",
    "  scale_x_continuous(limits=c(0,21), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,500),  expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "  \n",
    "\n",
    "histogram\n",
    "\n",
    "# ggsave('figures/break_length_histogram.pdf', plot=histogram, width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(paste('Mean break length: ', mean(df_break$dayDiffForward)))\n",
    "print(paste('StD break length: ', sd(df_break$dayDiffForward)))\n",
    "print(paste('Median break length: ', median(df_break$dayDiffForward)))\n",
    "print(paste('Min break length: ', min(df_break$dayDiffForward)))\n",
    "print(paste('Max break length: ', max(df_break$dayDiffForward)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_diff_distribution <- df_noClose[df_noClose['userDiffPrev'] != TRUE, ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(paste('Min', min(df_diff_distribution$hourDiffReverse)))\n",
    "print(paste('Max', max(df_diff_distribution$hourDiffReverse)))\n",
    "print(paste('Mean', mean(df_diff_distribution$hourDiffReverse)))\n",
    "print(paste('StD', sd(df_diff_distribution$hourDiffReverse)))\n",
    "print(paste('Median', median(df_diff_distribution$hourDiffReverse)))\n",
    "print(paste('Median minutes', median(df_diff_distribution$minuteDiffReverse)))\n",
    "print(paste('IQR', IQR(df_diff_distribution$hourDiffReverse)))\n",
    "print(paste('Q3', quantile(df_diff_distribution$hourDiffReverse, probs=0.75,na.rm = TRUE)))\n",
    "print(paste('Mean + stds', mean(df_diff_distribution$hourDiffReverse) + how_many_sds*sd(df_diff_distribution$hourDiffReverse)))\n",
    "\n",
    "mean_timeDiff_hours = mean(df_diff_distribution$hourDiffReverse)\n",
    "std_timeDiff_hours = sd(df_diff_distribution$hourDiffReverse)\n",
    "break_cutoff_timeDiff_hours = mean_timeDiff_hours + how_many_sds*std_timeDiff_hours\n",
    "\n",
    "# IQR=IQR(dv,na.rm = TRUE)\n",
    "#   Quant_25=quantile(dv,probs=0.25,na.rm = TRUE)\n",
    "#   Quant_75=quantile(dv,probs=0.75,na.rm = TRUE)\n",
    "#   upper=Quant_75+Tukey_crit*IQR\n",
    "#   lower=Quant_25-Tukey_crit*IQR\n",
    "#   outlier_Tukey=ifelse(dv>upper,1,ifelse(dv<lower,1,0))\n",
    "#   # print(outlier_Tukey)\n",
    "#   as.numeric(paste(outlier_Tukey))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# head(df_break)\n",
    "# Create the histogram using ggplot2\n",
    "histogram <- ggplot(data = df_diff_distribution, aes(x = minuteDiffReverse)) +\n",
    "  geom_histogram(binwidth = 10, fill = two_color[1], color = \"black\", alpha = 0.7) +\n",
    "  geom_vline(xintercept=24*60, color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=48*60, color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=mean(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=median(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=mean(df_diff_distribution$minuteDiffReverse) + 2*sd(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'solid') +\n",
    "  labs(x = \"Time Gap Between Uses (minutes)\", y = \"Frequency\") +\n",
    "  scale_x_continuous(limits=c(NA,3000), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,630000), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=20, r=25, b=0, l=0, unit='pt'))\n",
    "  \n",
    "\n",
    "histogram\n",
    "\n",
    "# ggsave('figures/usage_gap_histogram.pdf', plot=histogram, width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Create the histogram using ggplot2\n",
    "histogram <- ggplot(data = df_diff_distribution, aes(x = minuteDiffReverse)) +\n",
    "  geom_histogram(binwidth = 10, fill = two_color[1], color = \"black\", alpha = 0.7) +\n",
    "  geom_vline(xintercept=24*60, color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=mean(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=median(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'dashed') +\n",
    "  labs(x = \"Time Gap Between Uses (minutes)\", y = \"Frequency\") +\n",
    "  scale_x_continuous(limits=c(100,1000), expand=c(0,0), breaks = c(100,500,1000)) +\n",
    "  scale_y_continuous(limits=c(0,17000), expand=c(0,0), breaks = pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=25, b=0, l=0, unit='pt'))\n",
    "  \n",
    "\n",
    "histogram\n",
    "\n",
    "# ggsave('figures/usage_gap_histogram_zoom.pdf', plot=histogram, width=10, height=3.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Create the histogram using ggplot2\n",
    "histogram <- ggplot(data = df_diff_distribution, aes(x = minuteDiffReverse)) +\n",
    "  geom_histogram(binwidth = 10, fill = two_color[1], color = \"black\", alpha = 0.7) +\n",
    "  geom_vline(xintercept=24*60, color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=mean(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'dashed') +\n",
    "  geom_vline(xintercept=median(df_diff_distribution$minuteDiffReverse), color='black', linetype = 'dashed') +\n",
    "  labs(x = \"Time Gap Between Uses (minutes)\", y = \"Frequency\") +\n",
    "  scale_x_continuous(limits=c(NA,50), expand=c(0,0), breaks = c(100,500,1000)) +\n",
    "#   scale_y_continuous(limits=c(0,17000), expand=c(0,0), breaks = pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=25, b=0, l=0, unit='pt'))\n",
    "  \n",
    "\n",
    "histogram\n",
    "\n",
    "# ggsave('figures/usage_gap_histogram_zoom.pdf', plot=histogram, width=10, height=3.5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First Use, Leadup to Breaks and Returns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Everyone's first use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_first_use <- df_entries_per_user_hours[df_entries_per_user_hours$hoursSinceStart <= break_cutoff_timeDiff_hours,]\n",
    "df_first_use_result <- df_entries_per_user_hours_result[df_entries_per_user_hours_result$hoursSinceStart <= break_cutoff_timeDiff_hours,]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# counted_firsts <- df %>%\n",
    "#   group_by(hoursSinceStart, continuedToApp) %>%\n",
    "#   summarise(avg_entries_per_user = n() /n_distinct(userIndex))\n",
    "\n",
    "\n",
    "# counted_firsts <- drop_na(counted_firsts)\n",
    "df_first_use_result$continuedToApp <- factor(df_first_use_result$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "first_use_plot <- ggplot(df_first_use_result, aes(x = hoursSinceStart, y = interactions_per_user, fill = continuedToApp)) +\n",
    "  geom_bar(stat = \"identity\") +\n",
    "  labs(x = \"Time Since First Use (Hours)\", y = \"App Open Attempts per User\", fill = \"User Action\") +\n",
    "  theme(legend.position = 'none') +\n",
    "  scale_x_continuous(limits=c(NA, break_cutoff_timeDiff_hours), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,10), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=30), plot.margin = margin(t=10, r=10, b=0, l=5, unit='pt'))\n",
    "\n",
    "# pre_break_plot <- pre_break_plot + theme(legend.position = 'none')\n",
    "first_use_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# grouped_firsts <- counted_firsts %>%\n",
    "#   group_by(hoursSinceStart, continuedToApp) %>%\n",
    "#   summarize(\n",
    "#     average_user_count = mean(total_user_count),\n",
    "#     summed_normalized_count = sum(normalized_count)\n",
    "#   ) %>%\n",
    "#   mutate(\n",
    "#     percentage = summed_normalized_count / sum(summed_normalized_count) * 100,\n",
    "#     opens_per_user = sum(summed_normalized_count)\n",
    "#   )\n",
    "# df_first_use$continuedToApp <- factor(df_first_use$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "# grouped_firsts_effectiveness <- grouped_firsts[grouped_firsts$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "# grouped_firsts_effectiveness <- grouped_firsts_effectiveness[grouped_firsts_effectiveness$hoursSinceStart >= 0, ]\n",
    "\n",
    "secondary_scale <- 10\n",
    "\n",
    "# Create a dual-axis plot\n",
    "firsts_line <- ggplot(df_first_use, aes(x = hoursSinceStart)) +\n",
    "  geom_line(aes(y = interactions_per_user, color = \"A\"), size = 1) +\n",
    "  geom_line(data=df_first_use_result[df_first_use_result$continuedToApp == 'didNotOpenApp',], aes(y = percentage/secondary_scale, color = \"B\"), size = 1) +\n",
    "  scale_color_manual(values = c(two_color[2], two_color[1]), labels=c(\"App Open Attempts per User\",\"Dismissed Attempts (%)\")) +\n",
    "  labs(x = \"Time Since First Use (Hours)\", y = \"App Open Attempts per User\", color = \"\") +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=30),  legend.position = \"bottom\", plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt')) +\n",
    "  scale_y_continuous(limits=c(0,10), expand=c(0,0), breaks=pretty_breaks(),\n",
    "    sec.axis = sec_axis(~.*secondary_scale, name = \"Dismissed Attempts (%)\", labels=percent_format(scale = 1))) +\n",
    "  scale_x_continuous(limits=c(NA, break_cutoff_timeDiff_hours), expand=c(0,0), breaks=pretty_breaks())\n",
    "\n",
    "firsts_line\n",
    "\n",
    "\n",
    "\n",
    "# ggsave('figures/attempts_and_effectiveness.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Define the combined exponential-sinusoidal function to fit\n",
    "combined_model <- function(x, A, B, C, D, E, F, G) {\n",
    "  A * exp(-B * x) + C * cos(2*pi/24*D * x + E) + F * x + G\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_combo <- nls(interactions_per_user ~ combined_model(hoursSinceStart, A, B, C, D, E, F, G),\n",
    "        data=df_first_use,\n",
    "        start = list(A = 5, B = 0.9, C = 0.75, D = 1, E = 0, F=0, G=1),\n",
    "        control = nls.control(maxiter = 10000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_combo)\n",
    "# R2nls(fit_combo)\n",
    "\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "exp_model <- function(x, A, B, C, D) {\n",
    "  A * exp(-B * x) + C * x + D\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_exp <- nls(interactions_per_user ~ exp_model(hoursSinceStart, A, B, C, D),\n",
    "        data=df_first_use,\n",
    "        start = list(A = 5, B = .7, C = 0, D=5),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "sinusoidal <- function(x, A, B, C) {\n",
    "  A * cos(2*pi/24 * x + B) + C\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_sin <- nls(interactions_per_user ~ sinusoidal(hoursSinceStart, A, B, C),\n",
    "        data=df_first_use,\n",
    "        start = list(A = 3, B = 0, C = 1),\n",
    "        control = nls.control(maxiter = 5000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_sin)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(interactions_per_user ~ A*hoursSinceStart + B,\n",
    "        data=df_first_use,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_lin)\n",
    "\n",
    "\n",
    "\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = df_first_use$hoursSinceStart,\n",
    "                        y = df_first_use$interactions_per_user,\n",
    "                        ymin=df_first_use$interactions_per_user-df_first_use$std_error_interactions,\n",
    "                        ymax=df_first_use$interactions_per_user+df_first_use$std_error_interactions)\n",
    "                        \n",
    "line_x <- seq(min(df_first_use$hoursSinceStart),max(df_first_use$hoursSinceStart),.25)\n",
    "line_data <- data.frame(x = line_x,\n",
    "                        ymin=line_x,\n",
    "                        ymax=line_x,\n",
    "                        Fitted_combo = predict(fit_combo, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = line_x))\n",
    " )\n",
    "\n",
    "\n",
    "\n",
    " \n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[2], alpha = 0.5)+\n",
    "  geom_point(color = two_color[2], alpha=1) + # Data points\n",
    "  # geom_line(aes(y = Fitted_lin), color = two_color[1], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_sin), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_exp), color = two_color[1], size=1) + # Fitted curve\n",
    "  geom_line(data=line_data, aes(x=x, y = Fitted_combo), color = two_color[1], size=1) + # Fitted curve\n",
    "  annotate(\"text\", x=3, y=6, label=expression(R^2 * \"=.970\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  labs(x = \"Time Since First Use (Hours)\", y = \"App Open Attempts per User\") +\n",
    "  scale_x_continuous(limits=c(0,break_cutoff_timeDiff_hours), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,7), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg\n",
    "\n",
    "# ggsave('figures/first_fit.pdf', width=7, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "j=anova(fit_lin, fit_sin,fit_combo)\n",
    "j\n",
    "k=anova(fit_exp,fit_combo)\n",
    "k\n",
    "BIC(fit_lin, fit_sin, fit_exp,fit_combo)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "j[3,'Sum Sq']/sum(j['Res.Sum Sq'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_fit <- df_first_use_result[df_first_use_result$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "combined_model <- function(x, A, B, C, D, E, F, G) {\n",
    "  A * exp(-B * x) + C * cos(2*pi/24*D * x + E) + F * x + G\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_combo <- nls(percentage ~ combined_model(hoursSinceStart, A, B, C, D, E, F, G),\n",
    "        data=df_fit,\n",
    "        start = list(A = 5, B = 0.9, C = 0.75, D = 1, E = 0, F=0, G=1),\n",
    "        control = nls.control(maxiter = 10000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_combo)\n",
    "R2nls(fit_combo)\n",
    "\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "exp_model <- function(x, A, B, C, D) {\n",
    "  A * exp(-B * x) + C * x + D\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_exp <- nls(percentage ~ exp_model(hoursSinceStart, A, B, C, D),\n",
    "        data=df_fit,\n",
    "        start = list(A = 5, B = .7, C = 0, D=5),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "sinusoidal <- function(x, A, B, C) {\n",
    "  A * cos(2*pi/24 * x + B) + C\n",
    "}\n",
    "\n",
    "summary(fit_exp)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_sin <- nls(percentage ~ sinusoidal(hoursSinceStart, A, B, C),\n",
    "        data=df_fit,\n",
    "        start = list(A = 3, B = 0, C = 1),\n",
    "        control = nls.control(maxiter = 5000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_sin)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(percentage ~ A*hoursSinceStart + B,\n",
    "        data=df_fit,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_lin)\n",
    "\n",
    "\n",
    "\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = df_fit$hoursSinceStart,\n",
    "                        y = df_fit$percentage,\n",
    "                        ymin=df_fit$percentage-df_fit$std_error_percentage,\n",
    "                        ymax=df_fit$percentage+df_fit$std_error_percentage)\n",
    "line_x <- seq(min(df_fit$hoursSinceStart),max(df_fit$hoursSinceStart),.5)\n",
    "line_data <- data.frame(x = line_x,\n",
    "                        ymin=line_x,\n",
    "                        ymax=line_x,\n",
    "                        Fitted_combo = predict(fit_combo, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = line_x))\n",
    ")\n",
    " \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[1], alpha = 0.5)+\n",
    "  geom_point(color = two_color[1], alpha=1) + # Data points\n",
    "  # geom_line(aes(y = Fitted_lin), color = two_color[1], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_sin), color = two_color[2], size=1) + # Fitted curve\n",
    "  geom_line(data=line_data, aes(x=x, y = Fitted_exp), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_combo), color = two_color[1], size=1) + # Fitted curve\n",
    "  annotate(\"text\", x=3, y=100, label=expression(R^2 * \"=.807\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  labs(x = \"Time Since First Use (Hours)\", y = \"Dismissed Attempts (%)\") +\n",
    "  scale_x_continuous(limits=c(1,break_cutoff_timeDiff_hours), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,100), expand=c(0,0), breaks=c(0,10,20,30,40)) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg  \n",
    "\n",
    "# ggsave('figures/first_percentage_fit.pdf', width=7, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "j=anova(fit_lin,fit_sin,fit_exp,fit_combo)\n",
    "j\n",
    "BIC(fit_lin,fit_sin,fit_exp,fit_combo)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "eta_sq <- j[3,'Sum Sq']/sum(j['Res.Sum Sq'])\n",
    "eta_sq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Breaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_breaks<-read_csv(\"./data_breaks.csv\")\n",
    "colnames(df_breaks)<-make.names(colnames(df_breaks))\n",
    "# head(df_breaks, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_breaks <- df_breaks[df_breaks$hoursSinceStart >= -break_cutoff_timeDiff_hours,]\n",
    "number_users_breaks <- length(unique(df_breaks$userIndex))\n",
    "\n",
    "counted_breaks <- df_breaks %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize(interactions_per_user = n()/number_users_breaks)\n",
    "\n",
    "counted_breaks <- drop_na(counted_breaks)\n",
    "counted_breaks$continuedToApp <- factor(counted_breaks$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "pre_break_plot <- ggplot(counted_breaks, aes(x = time_to_break, y = interactions_per_user, fill = continuedToApp)) +\n",
    "  geom_bar(stat = \"identity\") +\n",
    "  labs(x = \"Time Before Break (Hours)\", y = \"App Open Attempts per User\", fill = \"User Action\") +\n",
    "  theme(legend.position = 'none') +\n",
    "  scale_x_continuous(limits=c(- break_cutoff_timeDiff_hours ,NA), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,10), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=30), plot.margin = margin(t=10, r=10, b=0, l=5, unit='pt'))\n",
    "\n",
    "# pre_break_plot <- pre_break_plot + theme(legend.position = 'none')\n",
    "pre_break_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_breaks <- counted_breaks %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize(\n",
    "    summed_normalized_count = sum(interactions_per_user)\n",
    "  ) %>%\n",
    "  mutate(\n",
    "    percentage = summed_normalized_count / sum(summed_normalized_count) * 100,\n",
    "    opens_per_user = sum(summed_normalized_count)\n",
    "  )\n",
    "\n",
    "grouped_breaks_effectiveness <- grouped_breaks[grouped_breaks$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "grouped_breaks_effectiveness <- grouped_breaks_effectiveness[grouped_breaks_effectiveness$time_to_break <= 0, ]\n",
    "\n",
    "secondary_scale <- 10\n",
    "\n",
    "# Create a dual-axis plot\n",
    "break_line <- ggplot(grouped_breaks_effectiveness, aes(x = time_to_break)) +\n",
    "  geom_line(aes(y = opens_per_user, color = \"A\"), size = 1) +\n",
    "  geom_line(aes(y = percentage/secondary_scale, color = \"B\"), size = 1) +\n",
    "  scale_color_manual(values = c(two_color[2], two_color[1]), labels=c(\"App Open Attempts per User\",\"Dismissed Attempts (%)\")) +\n",
    "  labs(x = \"Time Before Break (Hours)\", y = \"App Open Attempts per User\", color = \"\") +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=30),  legend.position = \"bottom\", plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt')) +\n",
    "  scale_y_continuous(limits=c(0,10), expand=c(0,0), breaks=pretty_breaks(),\n",
    "    sec.axis = sec_axis(~.*secondary_scale, name = \"Dismissed Attempts (%)\", labels=percent_format(scale = 1))) +\n",
    "  scale_x_continuous(limits=c(-break_cutoff_timeDiff_hours,NA), expand=c(0,0), breaks=pretty_breaks())\n",
    "\n",
    "break_line\n",
    "\n",
    "# ggsave('figures/attempts_and_effectiveness.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_breaks <- counted_breaks %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize(\n",
    "    summed_normalized_count = sum(interactions_per_user)\n",
    "  ) %>%\n",
    "  mutate(\n",
    "    percentage = summed_normalized_count / sum(summed_normalized_count) * 100,\n",
    "    opens_per_user = sum(summed_normalized_count)\n",
    "  )\n",
    "\n",
    "# grouped_breaks_effectiveness <- grouped_breaks[grouped_breaks$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "grouped_breaks_effectiveness <- grouped_breaks_effectiveness[grouped_breaks_effectiveness$time_to_break <= 0, ]\n",
    "\n",
    "# count how many breaks each user had\n",
    "df_breaks_per_user <- df_break %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(breaks_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_breaks_errors <- df_breaks %>%\n",
    "  group_by(userIndex, time_to_break, continuedToApp) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, df_breaks_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_breaks_errors$interactions_per_user <- counted_breaks_errors$interactions / counted_breaks_errors$breaks_per_user\n",
    "\n",
    "counted_breaks_errors <- counted_breaks_errors %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize( std_error_interactions = sd(interactions_per_user)/sqrt(length(interactions_per_user))\n",
    "            )\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, grouped_breaks, by=c('time_to_break','continuedToApp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_breaks <- counted_breaks %>%\n",
    "  group_by(time_to_break) %>%\n",
    "  summarize(\n",
    "    interactions_per_user = sum(interactions_per_user)\n",
    "  ) \n",
    "\n",
    "# grouped_breaks_effectiveness <- grouped_breaks[grouped_breaks$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "grouped_breaks_effectiveness <- grouped_breaks_effectiveness[grouped_breaks_effectiveness$time_to_break <= 0, ]\n",
    "\n",
    "# count how many breaks each user had\n",
    "df_breaks_per_user <- df_break %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(breaks_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_breaks_errors <- df_breaks %>%\n",
    "  group_by(userIndex, time_to_break) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, df_breaks_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_breaks_errors$interactions_per_user <- counted_breaks_errors$interactions / counted_breaks_errors$breaks_per_user\n",
    "\n",
    "counted_breaks_errors <- counted_breaks_errors %>%\n",
    "  group_by(time_to_break) %>%\n",
    "  summarize( std_error_interactions = sd(interactions_per_user)/sqrt(length(interactions_per_user))\n",
    "            )\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, grouped_breaks, by=c('time_to_break'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "counted_breaks_errors$hoursSinceStart <- counted_breaks_errors$time_to_break\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "combined_model <- function(x, A, B, C, D, E, F, G) {\n",
    "  A * exp(B * x) + C * cos(2*pi/24*D * x + E) + F * x + G\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_combo <- nls(interactions_per_user ~ combined_model(hoursSinceStart, A, B, C, D, E, F, G),\n",
    "        data=counted_breaks_errors,\n",
    "        start = list(A = 6, B = .06, C = 1, D = 1, E = 0, F=0, G=1),\n",
    "        control = nls.control(maxiter = 10000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "summary(fit_combo)\n",
    "R2nls(fit_combo)\n",
    "\n",
    "\n",
    "# # Define the combined exponential-sinusoidal function to fit\n",
    "# exp_model <- function(x, A, B, C, D, E) {\n",
    "#   A * exp(B * x + C) + D * x + E\n",
    "# }\n",
    "\n",
    "# # Fit the model to the data\n",
    "# fit_exp <- nls(interactions_per_user ~ exp_model(hoursSinceStart, A, B, C, D, E),\n",
    "#         data=counted_breaks_errors,\n",
    "#         start = list(A = 6, B = .06, C=1, D = 0, E=1),\n",
    "#         control = nls.control(maxiter = 10000)\n",
    "#         )\n",
    "\n",
    "# # Print the summary of the fit\n",
    "# summary(fit_exp)\n",
    "\n",
    "sinusoidal <- function(x, A, B, C) {\n",
    "  A * cos(2*pi/24 * x + B) + C\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_sin <- nls(interactions_per_user ~ sinusoidal(hoursSinceStart, A, B, C),\n",
    "        data=counted_breaks_errors,\n",
    "        start = list(A = 3, B = 0, C = 1),\n",
    "        control = nls.control(maxiter = 5000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_sin)\n",
    "# R2nls(fit_sin)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(interactions_per_user ~ A*hoursSinceStart + B,\n",
    "        data=counted_breaks_errors,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_lin)\n",
    "# R2nls(fit_lin)\n",
    "\n",
    "\n",
    "# # # Create a data frame for plotting\n",
    "# plot_data <- data.frame(x = df_second_use$hoursSinceStart,\n",
    "#                         y = df_second_use$interactions_per_user,\n",
    "#                         ymin=df_second_use$interactions_per_user-df_second_use$std_error_interactions,\n",
    "#                         ymax=df_second_use$interactions_per_user+df_second_use$std_error_interactions,\n",
    "#                         Fitted_combo = predict(fit_exp, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart)),\n",
    "#                         Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart)),\n",
    "#                         Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart)),\n",
    "#                         Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart))\n",
    "#  )\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = counted_breaks_errors$hoursSinceStart,\n",
    "                        y = counted_breaks_errors$interactions_per_user,\n",
    "                        ymin=counted_breaks_errors$interactions_per_user-counted_breaks_errors$std_error_interactions,\n",
    "                        ymax=counted_breaks_errors$interactions_per_user+counted_breaks_errors$std_error_interactions)\n",
    "\n",
    "line_x <- seq(min(counted_breaks_errors$hoursSinceStart),max(counted_breaks_errors$hoursSinceStart),.25)\n",
    "line_data <- data.frame(x = line_x,\n",
    "                        ymin=line_x,\n",
    "                        ymax=line_x,\n",
    "                        Fitted_combo = predict(fit_combo, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        # Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = line_x))\n",
    " )\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[2], alpha = 0.5)+\n",
    "  geom_point(color = two_color[2], alpha=1) + # Data points\n",
    "  # geom_line(aes(y = Fitted_lin), color = two_color[1], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_sin), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_exp), color = two_color[1], size=1) + # Fitted curve\n",
    "  geom_line(data=line_data, aes(x=x, y = Fitted_combo), color = two_color[1], size=1) + # Fitted curve\n",
    "  labs(x = \"Time Before Break (Hours)\", y = \"App Open Attempts per User\") +\n",
    "  annotate(\"text\", x=-45, y=6, label=expression(R^2 * \"=.978\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  scale_x_continuous(limits=c(-break_cutoff_timeDiff_hours,1), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,7), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg\n",
    "\n",
    "# ggsave('figures/break_fit.pdf', width=7, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "j=anova(fit_lin, fit_sin, fit_combo)\n",
    "j\n",
    "BIC(fit_lin, fit_sin, fit_combo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "eta_sq <- j[3,'Sum Sq']/sum(j['Res.Sum Sq'])\n",
    "eta_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_breaks <- counted_breaks %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize(\n",
    "    summed_normalized_count = sum(interactions_per_user)\n",
    "  ) %>%\n",
    "  mutate(\n",
    "    percentage = summed_normalized_count / sum(summed_normalized_count) * 100,\n",
    "    opens_per_user = sum(summed_normalized_count)\n",
    "  )\n",
    "\n",
    "# grouped_breaks_effectiveness <- grouped_breaks[grouped_breaks$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "grouped_breaks_effectiveness <- grouped_breaks_effectiveness[grouped_breaks_effectiveness$time_to_break <= 0, ]\n",
    "\n",
    "# count how many breaks each user had\n",
    "df_breaks_per_user <- df_break %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(breaks_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_breaks_errors <- df_breaks %>%\n",
    "  group_by(userIndex, time_to_break, continuedToApp) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, df_breaks_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_breaks_errors$interactions_per_user <- counted_breaks_errors$interactions / counted_breaks_errors$breaks_per_user\n",
    "\n",
    "counted_breaks_errors <- counted_breaks_errors %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize( std_error_interactions = sd(interactions_per_user)/sqrt(length(interactions_per_user))\n",
    "            )\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, grouped_breaks, by=c('time_to_break','continuedToApp'))\n",
    "counted_breaks_errors$squared_error_ratio = (counted_breaks_errors$std_error_interactions / counted_breaks_errors$opens_per_user)^2\n",
    "counted_breaks_errors <- counted_breaks_errors %>% \n",
    "  group_by(time_to_break) %>%\n",
    "  mutate(error_sum = sum(std_error_interactions),\n",
    "        interaction_sum = sum(opens_per_user)) %>%\n",
    "  ungroup()\n",
    "\n",
    "counted_breaks_errors$std_error_percentage <- counted_breaks_errors$percentage * sqrt(counted_breaks_errors$squared_error_ratio + (counted_breaks_errors$error_sum ^2) / (counted_breaks_errors$interaction_sum ^2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_fit <- counted_breaks_errors[counted_breaks_errors$continuedToApp == 'didNotOpenApp',]\n",
    "df_fit$hoursSinceStart <- df_fit$time_to_break\n",
    "\n",
    "# combined_model <- function(x, A, B, C, D, E, F, G) {\n",
    "#   A * exp(B * x) + C * cos(2*pi/24*D * x + E) + F * x + G\n",
    "# }\n",
    "\n",
    "# # Fit the model to the data\n",
    "# fit_combo <- nls(percentage ~ combined_model(hoursSinceStart, A, B, C, D, E, F, G),\n",
    "#         data=df_fit,\n",
    "#         start = list(A = 38, B = .02, C = 1, D = 1, E = 0, F=0, G=1),\n",
    "#         control = nls.control(maxiter = 10000)\n",
    "#         )\n",
    "\n",
    "# # Print the summary of the fit\n",
    "# summary(fit_combo)\n",
    "# # R2nls(fit_combo)\n",
    "\n",
    "\n",
    "# # # Define the combined exponential-sinusoidal function to fit\n",
    "# exp_model <- function(x, A, B, C, D, E) {\n",
    "#   A * exp(B * x + C) + D * x + E\n",
    "# }\n",
    "\n",
    "# # Fit the model to the data\n",
    "# fit_exp <- nls(percentage ~ exp_model(hoursSinceStart, A, B, C, D, E),\n",
    "#         data=df_fit,\n",
    "#         start = list(A = 38, B = .002, C=1, D = 0, E=1),\n",
    "#         control = nls.control(maxiter = 10000)\n",
    "#         )\n",
    "\n",
    "# # Print the summary of the fit\n",
    "# summary(fit_exp)\n",
    "\n",
    "sinusoidal <- function(x, A, B, C) {\n",
    "  A * cos(2*pi/24 * x + B) + C\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_sin <- nls(percentage ~ sinusoidal(hoursSinceStart, A, B, C),\n",
    "        data=df_fit,\n",
    "        start = list(A = 5, B = 0, C = 40),\n",
    "        control = nls.control(maxiter = 5000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_sin)\n",
    "# R2nls(fit_sin)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(percentage ~ A*hoursSinceStart + B,\n",
    "        data=df_fit,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_lin)\n",
    "# R2nls(fit_lin)\n",
    "\n",
    "\n",
    "\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = df_fit$hoursSinceStart,\n",
    "                        y = df_fit$percentage,\n",
    "                        ymin=df_fit$percentage-df_fit$std_error_percentage,\n",
    "                        ymax=df_fit$percentage+df_fit$std_error_percentage,\n",
    "                        # Fitted_combo = predict(fit_combo, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        # Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = df_fit$hoursSinceStart)),\n",
    "                        Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = df_fit$hoursSinceStart))\n",
    ")\n",
    " \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[1], alpha = 0.5)+\n",
    "  geom_point(color = two_color[1], alpha=1) + # Data points\n",
    "  # geom_smooth(method = \"lm\", se=FALSE, color=\"black\", formula = y ~ x) +\n",
    "  geom_line(aes(y = Fitted_lin), color = two_color[2], size=1, linetype=\"dashed\") + # Fitted curve\n",
    "  # geom_line(data=line_data, aes(x=x,y = Fitted_sin), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(data=line_data, aes(x=x, y = Fitted_exp), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_combo), color = two_color[1], size=1) + # Fitted curve\n",
    "  annotate(\"text\", x=-45, y=100, label=expression(R^2 * \"=.00647\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  labs(x = \"Time Before Break (Hours)\", y = \"Dismissed Attempts (%)\") +\n",
    "  scale_x_continuous(limits=c(-break_cutoff_timeDiff_hours,1), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,100), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg  \n",
    "\n",
    "# ggsave('figures/break_percentage_fit.pdf', width=7, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "anova(fit_sin, fit_lin )\n",
    "BIC(fit_lin, fit_sin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "lm_model <- lm(percentage ~ hoursSinceStart, data=df_fit)\n",
    "summary(lm_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "glmer_model <- lmer(percentage ~ hoursSinceStart + (1|summed_normalized_count), data=df_fit)\n",
    "summary(glmer_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(report(glmer_model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# count how many breaks each user had\n",
    "df_breaks_per_user <- df_break %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(breaks_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_breaks_errors <- df_breaks %>%\n",
    "  group_by(userIndex, time_to_break, continuedToApp) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_breaks_errors <- left_join(counted_breaks_errors, df_breaks_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_breaks_errors$interactions_per_user <- counted_breaks_errors$interactions / counted_breaks_errors$breaks_per_user\n",
    "\n",
    "counted_breaks_errors <- counted_breaks_errors %>%\n",
    "  group_by(time_to_break, continuedToApp) %>%\n",
    "  summarize(mean_interactions = mean(interactions_per_user),\n",
    "            std_error_interactions = sd(interactions_per_user)/sqrt(length(interactions_per_user))\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_returns<-read_csv(\"./data_returns.csv\")\n",
    "colnames(df_returns)<-make.names(colnames(df_returns))\n",
    "# head(df_returns, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_returns <- df_returns[df_returns$hoursSinceStart <= break_cutoff_timeDiff_hours,]\n",
    "number_users_returns <- length(unique(df_returns$userIndex))\n",
    "\n",
    "counted_returns <- df_returns %>%\n",
    "  group_by(time_after_return, continuedToApp) %>%\n",
    "  summarize(interactions_per_user = n()/number_users_returns)\n",
    "\n",
    "counted_returns$continuedToApp <- factor(counted_returns$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "post_return_plot <- ggplot(counted_returns, aes(x = time_after_return, y = interactions_per_user, fill = continuedToApp)) +\n",
    "  geom_bar(stat = \"identity\") +\n",
    "  labs(x = \"Hours After Return from Break\", y = \"App Open Attempts per User\", fill = \"User Action\") +\n",
    "  scale_x_continuous(limits=c(NA,break_cutoff_timeDiff_hours), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,10), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=30), plot.margin = margin(t=10, r=10, b=0, l=0, unit='pt'))\n",
    "\n",
    "# post_return_plot <- post_return_plot + theme(legend.position = 'none')\n",
    "post_return_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_returns <- counted_returns %>%\n",
    "  group_by(time_after_return, continuedToApp) %>%\n",
    "  summarize(\n",
    "    summed_normalized_count = sum(interactions_per_user)\n",
    "  ) %>%\n",
    "  mutate(\n",
    "    percentage = summed_normalized_count / sum(summed_normalized_count) * 100,\n",
    "    opens_per_user = sum(summed_normalized_count)\n",
    "  )\n",
    "\n",
    "grouped_returns_effectiveness <- grouped_returns[grouped_returns$continuedToApp == 'didNotOpenApp',]\n",
    "\n",
    "grouped_returns_effectiveness <- grouped_returns_effectiveness[grouped_returns_effectiveness$time_after_return >= 0, ]\n",
    "\n",
    "# Create a dual-axis plot\n",
    "return_line <- ggplot(grouped_returns_effectiveness, aes(x = time_after_return)) +\n",
    "  geom_line(aes(y = opens_per_user, color = \"A\"), size = 1) +\n",
    "  geom_line(aes(y = percentage/secondary_scale, color = \"B\"), size = 1) +\n",
    "  scale_color_manual(values = c(two_color[2], two_color[1]), labels=c(\"App Open Attempts per User\",\"Dismissed Attempts (%)\")) +\n",
    "  labs(x = \"Hours After Return from Break\", y = \"App Open Attempts per User\", color = \"\") +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=30),  legend.position = \"none\", plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt')) +\n",
    "  scale_y_continuous(limits=c(0,10), expand=c(0,0), breaks=pretty_breaks(),\n",
    "    sec.axis = sec_axis(~.*secondary_scale, name = \"Dismissed Attempts (%)\", labels = percent_format(scale = 1))) +\n",
    "  scale_x_continuous(limits=c(NA,break_cutoff_timeDiff_hours), expand=c(0,0), breaks=pretty_breaks())\n",
    "\n",
    "return_line\n",
    "\n",
    "# ggsave('figures/attempts_and_effectiveness.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_returns <- counted_returns %>%\n",
    "  group_by(time_after_return) %>%\n",
    "  summarize(\n",
    "    interactions_per_user = sum(interactions_per_user)\n",
    "  ) \n",
    "\n",
    "\n",
    "# count how many breaks each user had\n",
    "df_returns_per_user <- df_return %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(returns_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_returns_errors <- df_returns %>%\n",
    "  group_by(userIndex, time_after_return) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_returns_errors <- left_join(counted_returns_errors, df_returns_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_returns_errors$interactions_per_user <- counted_returns_errors$interactions / counted_returns_errors$returns_per_user\n",
    "\n",
    "counted_returns_errors <- counted_returns_errors %>%\n",
    "  group_by(time_after_return) %>%\n",
    "  summarize( std_error_interactions = sd(interactions_per_user)/sqrt(length(interactions_per_user))\n",
    "            )\n",
    "counted_returns_errors <- left_join(counted_returns_errors, grouped_returns, by=c('time_after_return'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "counted_returns_errors$hoursSinceStart <- counted_returns_errors$time_after_return\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "combined_model <- function(x, A, B, C, D, E, F, G) {\n",
    "  A * exp(-B * x) + C * cos(2*pi/24*D * x + E) + F * x + G\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_combo <- nls(interactions_per_user ~ combined_model(hoursSinceStart, A, B, C, D, E, F, G),\n",
    "        data=counted_returns_errors,\n",
    "        start = list(A = 1, B = 1, C = 1, D = 1, E = 0, F=0, G=1),\n",
    "        control = nls.control(maxiter = 10000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_combo)\n",
    "R2nls(fit_combo)\n",
    "\n",
    "\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "# exp_model <- function(x, A, B, C, D, E) {\n",
    "#   A * exp(-B * x - C) + D * x + E\n",
    "# }\n",
    "\n",
    "# # Fit the model to the data\n",
    "# fit_exp <- nls(interactions_per_user ~ exp_model(hoursSinceStart, A, B, C, D, E),\n",
    "#         data=counted_returns_errors,\n",
    "#         start = list(A = 1.1, B = 20, C=48, D = 0, E=1.1),\n",
    "#         control = nls.control(maxiter = 10000)\n",
    "#         )\n",
    "\n",
    "# # Print the summary of the fit\n",
    "# summary(fit_exp)\n",
    "\n",
    "sinusoidal <- function(x, A, B, C) {\n",
    "  A * cos(2*pi/24 * x + B) + C\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_sin <- nls(interactions_per_user ~ sinusoidal(hoursSinceStart, A, B, C),\n",
    "        data=counted_returns_errors,\n",
    "        start = list(A = 3, B = 0, C = 1),\n",
    "        control = nls.control(maxiter = 5000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_sin)\n",
    "# R2nls(fit_sin)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(interactions_per_user ~ A*hoursSinceStart + B,\n",
    "        data=counted_returns_errors,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "fit_exp <- nls(interactions_per_user ~ A*hoursSinceStart + B,\n",
    "        data=counted_returns_errors,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_lin)\n",
    "# R2nls(fit_lin)\n",
    "\n",
    "\n",
    "# # # Create a data frame for plotting\n",
    "# plot_data <- data.frame(x = df_second_use$hoursSinceStart,\n",
    "#                         y = df_second_use$interactions_per_user,\n",
    "#                         ymin=df_second_use$interactions_per_user-df_second_use$std_error_interactions,\n",
    "#                         ymax=df_second_use$interactions_per_user+df_second_use$std_error_interactions,\n",
    "#                         Fitted_combo = predict(fit_exp, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart)),\n",
    "#                         Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart)),\n",
    "#                         Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart)),\n",
    "#                         Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = df_second_use$hoursSinceStart))\n",
    "#  )\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = counted_returns_errors$hoursSinceStart,\n",
    "                        y = counted_returns_errors$interactions_per_user,\n",
    "                        ymin=counted_returns_errors$interactions_per_user-counted_returns_errors$std_error_interactions,\n",
    "                        ymax=counted_returns_errors$interactions_per_user+counted_returns_errors$std_error_interactions)\n",
    "\n",
    "line_x <- seq(min(counted_returns_errors$hoursSinceStart),max(counted_returns_errors$hoursSinceStart),.25)\n",
    "line_data <- data.frame(x = line_x,\n",
    "                        ymin=line_x,\n",
    "                        ymax=line_x,\n",
    "                        Fitted_combo = predict(fit_combo, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        # Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = line_x))\n",
    " )\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[2], alpha = 0.5)+\n",
    "  geom_point(color = two_color[2], alpha=1) + # Data points\n",
    "  # geom_line(aes(y = Fitted_lin), color = two_color[1], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_sin), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(aes(y = Fitted_exp), color = two_color[1], size=1) + # Fitted curve\n",
    "  geom_line(data=line_data, aes(x=x, y = Fitted_combo), color = two_color[1], size=1) + # Fitted curve\n",
    "  labs(x = \"Time After Return From Break (Hours)\", y = \"App Open Attempts per User\") +\n",
    "  annotate(\"text\", x=3, y=6, label=expression(R^2 * \"=.912\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  scale_x_continuous(limits=c(0,break_cutoff_timeDiff_hours), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,7), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg\n",
    "\n",
    "# ggsave('figures/return_fit.pdf', width=7, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "j=anova(fit_lin, fit_sin, fit_combo)\n",
    "j\n",
    "AIC(fit_lin, fit_sin,  fit_combo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "eta_sq <- j[3,'Sum Sq']/sum(j['Res.Sum Sq'])\n",
    "eta_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_returns <- counted_returns %>%\n",
    "  group_by(time_after_return, continuedToApp) %>%\n",
    "  summarize(\n",
    "    summed_normalized_count = sum(interactions_per_user)\n",
    "  ) %>%\n",
    "  mutate(\n",
    "    percentage = summed_normalized_count / sum(summed_normalized_count) * 100,\n",
    "    opens_per_user = sum(summed_normalized_count)\n",
    "  )\n",
    "\n",
    "# count how many breaks each user had\n",
    "df_returns_per_user <- df_return %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(returns_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_returns_errors <- df_returns %>%\n",
    "  group_by(userIndex, time_after_return, continuedToApp) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_returns_errors <- left_join(counted_returns_errors, df_returns_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_returns_errors$interactions_per_user <- counted_returns_errors$interactions / counted_returns_errors$returns_per_user\n",
    "\n",
    "counted_returns_errors <- counted_returns_errors %>%\n",
    "  group_by(time_after_return, continuedToApp) %>%\n",
    "  summarize( std_error_interactions = sd(interactions_per_user)/sqrt(length(interactions_per_user))\n",
    "            )\n",
    "counted_returns_errors <- left_join(counted_returns_errors, grouped_returns, by=c('time_after_return','continuedToApp'))\n",
    "counted_returns_errors$squared_error_ratio = (counted_returns_errors$std_error_interactions / counted_returns_errors$opens_per_user)^2\n",
    "counted_returns_errors <- counted_returns_errors %>% \n",
    "  group_by(time_after_return) %>%\n",
    "  mutate(error_sum = sum(std_error_interactions),\n",
    "        interaction_sum = sum(opens_per_user)) %>%\n",
    "  ungroup()\n",
    "\n",
    "counted_returns_errors$std_error_percentage <- counted_returns_errors$percentage * sqrt(counted_returns_errors$squared_error_ratio + (counted_returns_errors$error_sum ^2) / (counted_returns_errors$interaction_sum ^2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_fit <- counted_returns_errors[counted_returns_errors$continuedToApp == 'didNotOpenApp',]\n",
    "df_fit$hoursSinceStart <- df_fit$time_after_return\n",
    "\n",
    "## Define the combined exponential-sinusoidal function to fit\n",
    "combined_model <- function(x, A, B, C, D, E, F, G) {\n",
    "  A * exp(-B * x) + C * cos(2*pi/24*D * x + E) + F * x + G\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_combo <- nls(percentage ~ combined_model(hoursSinceStart, A, B, C, D, E, F, G),\n",
    "        data=df_fit,\n",
    "        start = list(A = 40, B = 0.9, C = 0.75, D = 1, E = 0, F=0, G=1),\n",
    "        control = nls.control(maxiter = 10000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_combo)\n",
    "# R2nls(fit_combo)\n",
    "\n",
    "# Define the combined exponential-sinusoidal function to fit\n",
    "exp_model <- function(x, A, B, C, D) {\n",
    "  A * exp(-B * x) + C * x + D\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_exp <- nls(percentage ~ exp_model(hoursSinceStart, A, B, C, D),\n",
    "        data=df_fit,\n",
    "        start = list(A = 17, B = 0.1, C = 0, D=48),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "\n",
    "# summary(fit_exp)\n",
    "R2nls(fit_exp)\n",
    "\n",
    "sinusoidal <- function(x, A, B, C) {\n",
    "  A * cos(2*pi/24 * x + B) + C\n",
    "}\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_sin <- nls(percentage ~ sinusoidal(hoursSinceStart, A, B, C),\n",
    "        data=df_fit,\n",
    "        start = list(A = 3, B = 0, C = 1),\n",
    "        control = nls.control(maxiter = 5000)\n",
    "        )\n",
    "\n",
    "# Print the summary of the fit\n",
    "# summary(fit_sin)\n",
    "# R2nls(fit_sin)\n",
    "\n",
    "# Fit the model to the data\n",
    "fit_lin <- nls(percentage ~ A*hoursSinceStart + B,\n",
    "        data=df_fit,\n",
    "        start = list(A = -1, B = 16),\n",
    "        control = nls.control(maxiter = 1000)\n",
    "        )\n",
    "# R2nls(fit_lin)\n",
    "\n",
    "\n",
    "# # Create a data frame for plotting\n",
    "plot_data <- data.frame(x = df_fit$hoursSinceStart,\n",
    "                        y = df_fit$percentage,\n",
    "                        ymin=df_fit$percentage-df_fit$std_error_percentage,\n",
    "                        ymax=df_fit$percentage+df_fit$std_error_percentage)\n",
    "line_x <- seq(min(df_fit$hoursSinceStart),max(df_fit$hoursSinceStart),.5)\n",
    "line_data <- data.frame(x = line_x,\n",
    "                        ymin=line_x,\n",
    "                        ymax=line_x,\n",
    "                        Fitted_combo = predict(fit_combo, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_exp = predict(fit_exp, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_sin = predict(fit_sin, newdata = data.frame(hoursSinceStart = line_x)),\n",
    "                        Fitted_lin = predict(fit_lin, newdata = data.frame(hoursSinceStart = line_x))\n",
    ")\n",
    " \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Create a ggplot\n",
    "gg <- ggplot(plot_data, aes(x = x, y = y, ymin=ymin, ymax=ymax)) +\n",
    "  geom_errorbar(width=0.2, color=two_color[1], alpha = 0.5)+\n",
    "  geom_point(color = two_color[1], alpha=1) + # Data points\n",
    "  # geom_line(data=line_data, aes(x=x,y = Fitted_lin), color = two_color[1], size=1) + # Fitted curve\n",
    "  # geom_line(data=line_data, aes(x=x,y = Fitted_sin), color = two_color[2], size=1) + # Fitted curve\n",
    "  geom_line(data=line_data, aes(x=x, y = Fitted_exp), color = two_color[2], size=1) + # Fitted curve\n",
    "  # geom_line(data=line_data, aes(x=x, y = Fitted_combo), color = two_color[2], size=1) + # Fitted curve\n",
    "  annotate(\"text\", x=3, y=100, label=expression(R^2 * \"=.370\"), parse=TRUE, vjust=1, hjust=0, size=7) +\n",
    "  labs(x = \"Time After Return From Break (Hours)\", y = \"Dismissed Attempts (%)\") +\n",
    "  scale_x_continuous(limits=c(0,break_cutoff_timeDiff_hours), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,100), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=12, b=0, l=0, unit='pt'))\n",
    "\n",
    "# Print the ggplot\n",
    "gg  \n",
    "\n",
    "# ggsave('figures/return_percentage_fit.pdf', width=7, height=7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "j=anova(fit_lin, fit_sin, fit_exp, fit_combo)\n",
    "j\n",
    "BIC(fit_lin, fit_sin, fit_exp, fit_combo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "eta_sq <- j[3,'Sum Sq']/sum(j['Res.Sum Sq'])\n",
    "eta_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# grid_plot <- grid.arrange(pre_break_plot, post_return_plot, pre_dropout_plot, ncol = 3)\n",
    "# final_plot <- arrangeGrob(grid_plot, bottom = legend_table)\n",
    "# grid.draw(final_plot)\n",
    "grid_plot <- grid.arrange(first_use_plot + theme(legend.position = \"none\"),\n",
    "                            pre_break_plot + theme(legend.position = \"none\", axis.title.y = element_blank()),\n",
    "                            post_return_plot + theme(legend.position = \"none\", axis.title.y = element_blank()), ncol = 3)\n",
    "\n",
    "temp <- pre_break_plot + theme(legend.position = \"bottom\", legend.direction = \"horizontal\",text = element_text(size=25))\n",
    "# Create a custom legend\n",
    "legend_grid <- get_legend(temp)\n",
    "\n",
    "# Add the legend to the combined plot\n",
    "# final_plot <- arranged_plots + theme(legend.position = \"none\")\n",
    "final_plot <- arrangeGrob(grid_plot, legend_grid, nrow=2, heights = c(0.9, 0.1))\n",
    "\n",
    "final_plot\n",
    "\n",
    "# ggsave('figures/grid_break_return_dropout.pdf', plot=final_plot, width=25, height=8, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# grid_plot <- grid.arrange(pre_break_plot, post_return_plot, pre_dropout_plot, ncol = 3)\n",
    "# final_plot <- arrangeGrob(grid_plot, bottom = legend_table)\n",
    "# grid.draw(final_plot)\n",
    "grid_plot <- grid.arrange(firsts_line + theme(legend.position = \"none\", axis.title.y.right = element_blank()),\n",
    "                            break_line + theme(legend.position = \"none\", axis.title.y = element_blank()),\n",
    "                            return_line + theme(legend.position = \"none\", axis.title.y = element_blank()), ncol = 3)\n",
    "\n",
    "temp <- break_line + theme(legend.position = \"bottom\", legend.direction = \"horizontal\", text = element_text(size=25))\n",
    "# Create a custom legend\n",
    "legend_grid <- get_legend(temp)\n",
    "\n",
    "# Add the legend to the combined plot\n",
    "# final_plot <- arranged_plots + theme(legend.position = \"none\")\n",
    "final_plot <- arrangeGrob(grid_plot, legend_grid, nrow=2, heights = c(0.9, 0.1))\n",
    "\n",
    "final_plot\n",
    "\n",
    "# ggsave('figures/grid_line_break_return_dropout.pdf', plot=final_plot, width=25, height=8, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_start <- df_noClose[df_noClose$hoursSinceStart <= break_cutoff_timeDiff_hours,]\n",
    "counted_firsts_box <- df_start %>%\n",
    "  group_by(userIndex, continuedToApp) %>%\n",
    "  summarize(interactions_per_user = n())\n",
    "\n",
    "# count how many breaks each user had\n",
    "df_breaks_per_user <- df_break %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(breaks_per_user = n())\n",
    "\n",
    "# count total interactions per user\n",
    "counted_breaks_box <- df_breaks %>%\n",
    "  group_by(userIndex, continuedToApp) %>%\n",
    "  summarize(interactions = n())\n",
    "\n",
    "# join the total interactions and number of breaks\n",
    "counted_breaks_box <- left_join(counted_breaks_box, df_breaks_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_breaks_box$interactions_per_user <- counted_breaks_box$interactions / counted_breaks_box$breaks_per_user\n",
    "\n",
    "\n",
    "# count how many returns each user had\n",
    "df_returns_per_user <- df_return %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarize(returns_per_user = n())\n",
    "\n",
    "counted_returns_box <- df_returns %>%\n",
    "  group_by(userIndex, continuedToApp) %>%\n",
    "  summarize(interactions = n()) \n",
    "\n",
    "# join the total interactions and number of returns\n",
    "counted_returns_box <- left_join(counted_returns_box, df_returns_per_user, by='userIndex')\n",
    "# create a new column where the total interactions is normalized by the number of breaks, so that it represents their use per break\n",
    "counted_returns_box$interactions_per_user <- counted_returns_box$interactions / counted_returns_box$returns_per_user\n",
    "\n",
    "counted_firsts_box$identifier <- as.factor(1) #\"First\"\n",
    "counted_breaks_box$identifier <- as.factor(2) #\"Break\"\n",
    "counted_returns_box$identifier <- as.factor(3) #\"Return\"\n",
    "\n",
    "\n",
    "counted_firsts_box$identifier_text <- as.factor('First') #\"First\"\n",
    "counted_breaks_box$identifier_text <- as.factor('Break') #\"Break\"\n",
    "counted_returns_box$identifier_text <- as.factor('Return') #\"Return\"\n",
    "\n",
    "counted_combined <- rbind(counted_firsts_box[c('userIndex','continuedToApp', 'interactions_per_user', 'identifier', 'identifier_text')], counted_breaks_box[c('userIndex','continuedToApp', 'interactions_per_user', 'identifier', 'identifier_text')], counted_returns_box[c('userIndex','continuedToApp', 'interactions_per_user', 'identifier', 'identifier_text')])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "counted_combined$continuedToApp <- factor(counted_combined$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "combined_box <- ggplot(data=counted_combined, aes(x=identifier_text, y = interactions_per_user, fill = continuedToApp)) +\n",
    "  geom_boxplot() +\n",
    "  labs(y = \"App Open Attempts per User\", fill = \"User Action\") +\n",
    "  # scale_x_discrete(labels=c(\"After First Use\",\"Before Break\", \"After Return\"), ) +\n",
    "  # scale_y_continuous(limits=c(0,40), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = 'bottom', legend.direction='horizontal',axis.title.x=element_blank(), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'))\n",
    "\n",
    "combined_box\n",
    "\n",
    "# ggsave('figures/user_action_box.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "counted_combined$continuedToApp <- factor(counted_combined$continuedToApp, levels=c('openedApp', 'didNotOpenApp'))\n",
    "\n",
    "combined_box <- ggplot(data=counted_combined, aes(x=identifier_text, y = interactions_per_user)) +\n",
    "  geom_boxplot(fill=plot_colors[1:3]) +\n",
    "  labs(y = \"App Open Attempts per User\") +\n",
    "  # scale_x_discrete(labels=c(\"After First Use\",\"Before Break\", \"After Return\"), ) +\n",
    "  scale_y_continuous(limits=c(0,40), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "#   scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  # scale_color_manual(values = two_color) + \n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = 'bottom', legend.direction='horizontal',axis.title.x=element_blank(), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'))\n",
    "\n",
    "combined_box\n",
    "\n",
    "# ggsave('figures/user_action_box.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "glmer_model <- lmer(interactions_per_user ~ identifier_text + (1 | userIndex), data=counted_combined)\n",
    "summary(glmer_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(report(glmer_model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "counted_combined_box <- counted_combined\n",
    "outliers_combined <- tukey_detect(dv=counted_combined_box$interactions_per_user)\n",
    "counted_combined_box[outliers_combined ==1,]$interactions_per_user <- as.integer(mean(counted_combined_box$interactions_per_user))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "combined_box <- ggplot(data=counted_combined_box, aes(x=identifier_text, y = interactions_per_user, fill = continuedToApp)) +\n",
    "  geom_boxplot() +\n",
    "  labs(y = \"App Open Attempts per User\", fill = \"User Action\") +\n",
    "  scale_x_discrete(labels=c(\"First Use\",\"Pre-Break\", \"Post-Return\"), ) +\n",
    "  # scale_y_continuous(limits=c(0,1), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  scale_fill_discrete(labels=c('Continued to App', 'Did not Continue')) +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = 'bottom', legend.direction='horizontal',axis.title.x=element_blank(), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'))\n",
    "\n",
    "combined_box\n",
    "\n",
    "# ggsave('figures/user_action_box.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "glmer_model <- lmer(interactions_per_user ~ identifier_text*continuedToApp + (1 | userIndex), data=counted_combined_box)\n",
    "summary(glmer_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(report(glmer_model))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical Testing for Pre-Post Break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "counted_breaks$identifier <- 'break'\n",
    "counted_returns$identifier <- 'return'\n",
    "temp <- rbind(counted_breaks, counted_returns)\n",
    "\n",
    "m <- aov(interactions_per_user ~ identifier * continuedToApp, data = temp)\n",
    "summary(m) # sig\n",
    "TukeyHSD(m) # sig Break-return and dropout-return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "temp$identifier <- as.factor(temp$identifier)\n",
    "m = art(interactions_per_user ~ identifier*continuedToApp, data = temp)\n",
    "anova(m) # sig\n",
    "art.con(m, ~identifier*continuedToApp) # sig Break-return and dropout-return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dickey Fuller Test for Stationary Time Series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "grouped_breaks_effectiveness$condition <- \"Break\"\n",
    "grouped_breaks_effectiveness$time <- grouped_breaks_effectiveness$time_to_break\n",
    "grouped_breaks_effectiveness <- grouped_breaks_effectiveness[!is.na(grouped_breaks_effectiveness$opens_per_user),]\n",
    "grouped_returns_effectiveness$condition <- \"Return\"\n",
    "grouped_returns_effectiveness$time <- grouped_returns_effectiveness$time_after_return\n",
    "grouped_returns_effectiveness <- grouped_returns_effectiveness[!is.na(grouped_returns_effectiveness$opens_per_user),]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(adf.test(df_first_use$interactions_per_user))\n",
    "print(adf.test(grouped_breaks_effectiveness$opens_per_user))\n",
    "print(adf.test(grouped_returns_effectiveness$opens_per_user))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(adf.test(df_first_use_result[df_first_use_result$continuedToApp == 'didNotOpenApp',]$percentage))\n",
    "print(adf.test(grouped_breaks_effectiveness$percentage))\n",
    "print(adf.test(grouped_returns_effectiveness$percentage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mann-Kendall Trend Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(MannKendall(df_first_use$interactions_per_user))\n",
    "print(MannKendall(grouped_breaks_effectiveness$opens_per_user))\n",
    "print(MannKendall(grouped_returns_effectiveness$opens_per_user))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(MannKendall(df_first_use_result[df_first_use_result$continuedToApp == 'didNotOpenApp',]$percentage))\n",
    "print(MannKendall(grouped_breaks_effectiveness$percentage))\n",
    "print(MannKendall(grouped_returns_effectiveness$percentage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pearson Correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "pearsonTest(df_first_use$hoursSinceStart, df_first_use$interactions_per_user)\n",
    "pearsonTest(grouped_breaks_effectiveness$time, grouped_breaks_effectiveness$opens_per_user)\n",
    "pearsonTest(grouped_returns_effectiveness$time, grouped_returns_effectiveness$opens_per_user)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Days of the Weeks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Create a histogram\n",
    "weekday_histogram <- df_noClose %>%\n",
    "  group_by(day_of_week) %>%\n",
    "  summarize(count = n())\n",
    "\n",
    "# Order days of the week properly\n",
    "days_ordered <- c(\"Monday\", \"Tuesday\", \"Wednesday\", \"Thursday\", \"Friday\", \"Saturday\",\"Sunday\")\n",
    "weekday_histogram$day_of_week <- factor(weekday_histogram$day_of_week, levels = days_ordered)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_combo <- df_noClose[df_noClose$lifeStatus %in% c(\"Break\",\"Return\"), ]\n",
    "\n",
    "df_combo$lifeStatus <- factor(df_combo$lifeStatus)\n",
    "\n",
    "\n",
    "weekday_histogram_combo <- df_combo %>%\n",
    "  group_by(day_of_week, lifeStatus) %>%\n",
    "  summarize(count = n())\n",
    "\n",
    "weekday_histogram_combo$day_of_week <- factor(weekday_histogram_combo$day_of_week, levels = days_ordered)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_alive <- df_noClose[df_noClose$lifeStatus %in% c(\"Alive\"), ]\n",
    "\n",
    "df_alive$lifeStatus <- \"Normal Use\"\n",
    "df_alive$lifeStatus <- factor(df_alive$lifeStatus)\n",
    "\n",
    "\n",
    "weekday_histogram_combo <- df_alive %>%\n",
    "  group_by(day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n()/length(unique(df_alive$userIndex)))\n",
    "\n",
    "weekday_histogram_combo$day_of_week <- factor(weekday_histogram_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Create the histogram plot\n",
    "ggplot(weekday_histogram_combo, aes(x = day_of_week)) +\n",
    "  geom_bar(stat = \"identity\", position=\"dodge\", aes(y = count_per_user, fill=lifeStatus)) +\n",
    "  labs(x = \"Day of the Week\", y = \"Count per User\", fill=\"Status\") +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_y_continuous(expand=c(0,0))+\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\")\n",
    "\n",
    "\n",
    "# ggsave('figures/day_of_week_alive.pdf', width=12, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_day_break <- df_noClose[df_noClose$lifeStatus %in% c(\"Alive\"), ]\n",
    "\n",
    "df_alive$lifeStatus <- \"Normal Use\"\n",
    "df_alive$lifeStatus <- factor(df_alive$lifeStatus)\n",
    "\n",
    "weekday_histogram_combo <- df_alive %>%\n",
    "  group_by(userIndex, day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n())\n",
    "\n",
    "df_max_weeks <- df %>%\n",
    "  group_by(userIndex) %>%\n",
    "  summarise(max_weeks = max(weeksSinceStart)) \n",
    "\n",
    "df_weekday <- left_join(weekday_histogram_combo, df_max_weeks, by='userIndex')\n",
    "df_weekday$normalized_count <- df_weekday$count_per_user / df_weekday$max_weeks\n",
    "\n",
    "\n",
    "\n",
    "df_weekday$day_of_week <- factor(df_weekday$day_of_week, levels = days_ordered)\n",
    "\n",
    "df_weekday_outlier <- df_weekday\n",
    "outliers_days <- tukey_detect(dv=df_weekday_outlier$normalized_count)\n",
    "df_weekday_outlier[outliers_days == 1,]$normalized_count <- mean(df_weekday_outlier[df_weekday_outlier$normalized_count != 'Inf',]$normalized_count)\n",
    "\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Create the histogram plot\n",
    "ggplot(df_weekday_outlier, aes(x = day_of_week, y = normalized_count, fill=lifeStatus)) +\n",
    "  geom_boxplot() +\n",
    "  stat_summary(fun = mean, geom = \"point\", shape = 4, size = 3, color = \"black\") +\n",
    "  labs(x = \"Day of the Week\", y = \"Count per User\", fill=\"Status\") +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_y_continuous(limits=c(0,50))+\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\")\n",
    "\n",
    "\n",
    "# ggsave('figures/day_of_week_box.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_day_break <- df_noClose[df_noClose$lifeStatus %in% c(\"Break\"), ]\n",
    "\n",
    "df_day_break$lifeStatus <- factor(df_day_break$lifeStatus)\n",
    "\n",
    "df_day_break_histogram_combo <- df_day_break %>%\n",
    "  group_by(userIndex, day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n()) %>%\n",
    "  group_by(day_of_week, lifeStatus) %>%\n",
    "  summarize(mean_interactions = mean(count_per_user),\n",
    "            std_error_interactions = sd(count_per_user)/sqrt(length(count_per_user))\n",
    "            )\n",
    "\n",
    "df_day_break_histogram_combo$day_of_week <- factor(df_day_break_histogram_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "df_day_return <- df_noClose[df_noClose$lifeStatus %in% c(\"Return\"), ]\n",
    "\n",
    "df_day_return$lifeStatus <- factor(df_day_return$lifeStatus)\n",
    "\n",
    "df_day_return_histogram_combo <- df_day_return %>%\n",
    "  group_by(userIndex, day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n()) %>%\n",
    "  group_by(day_of_week, lifeStatus) %>%\n",
    "  summarize(mean_interactions = mean(count_per_user),\n",
    "            std_error_interactions = sd(count_per_user)/sqrt(length(count_per_user))\n",
    "            )\n",
    "\n",
    "df_day_return_histogram_combo$day_of_week <- factor(df_day_return_histogram_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "df_day_break_return_histogram_combo <- rbind(df_day_break_histogram_combo, df_day_return_histogram_combo)\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Create the histogram plot\n",
    "ggplot(df_day_break_return_histogram_combo, aes(x = day_of_week,y = mean_interactions, fill=lifeStatus, ymin=mean_interactions-std_error_interactions, ymax=mean_interactions+std_error_interactions)) +\n",
    "  geom_bar( stat = \"identity\", position=\"dodge\", color='black') +\n",
    "  geom_errorbar(width=0.2, position=position_dodge(.9), color='black')+\n",
    "  labs(x = \"Day of the Week\", y = \"Normapized Count per User\", fill=\"Status\") +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_y_continuous(expand=c(0,0))+\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\")\n",
    "\n",
    "\n",
    "# ggsave('figures/day_of_week.pdf', width=20, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_day_break <- df_noClose[df_noClose$lifeStatus %in% c(\"Break\"), ]\n",
    "\n",
    "df_day_break$lifeStatus <- factor(df_day_break$lifeStatus)\n",
    "\n",
    "df_day_break_stats_combo <- df_day_break %>%\n",
    "  group_by(userIndex, day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n())\n",
    "\n",
    "# days_ordered <- c(\"Sunday\", \"Monday\", \"Tuesday\", \"Wednesday\", \"Thursday\", \"Friday\", \"Saturday\")\n",
    "days_ordered <- c(\"Wednesday\", \"Thursday\", \"Friday\", \"Saturday\",\"Sunday\", \"Monday\", \"Tuesday\")\n",
    "df_day_break_stats_combo$day_of_week <- factor(df_day_break_stats_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "\n",
    "df_day_return <- df_noClose[df_noClose$lifeStatus %in% c(\"Return\"), ]\n",
    "\n",
    "df_day_return$lifeStatus <- factor(df_day_return$lifeStatus)\n",
    "\n",
    "df_day_return_stats_combo <- df_day_return %>%\n",
    "  group_by(userIndex, day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n())\n",
    "\n",
    "df_day_return_stats_combo$day_of_week <- factor(df_day_return_stats_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "df_day_break_return_stats_combo <- rbind(df_day_break_stats_combo, df_day_return_stats_combo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "head(df_day_break_return_stats_combo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "glmer_model <- lmer(count_per_user ~ day_of_week + (1 | userIndex), data = df_day_break_stats_combo)\n",
    "summary(glmer_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(report(glmer_model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "glmer_model <- lmer(count_per_user ~ day_of_week + (1 | userIndex), data=df_day_return_stats_combo)\n",
    "summary(glmer_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(report(glmer_model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df_day_break <- df_noClose[df_noClose$lifeStatus %in% c(\"Break\"), ]\n",
    "\n",
    "df_day_break$lifeStatus <- factor(df_day_break$lifeStatus)\n",
    "\n",
    "\n",
    "df_day_break_histogram_combo <- df_day_break %>%\n",
    "  group_by(userIndex,day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n())\n",
    "\n",
    "df_day_break_histogram_combo$day_of_week <- factor(df_day_break_histogram_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "\n",
    "df_day_return <- df_noClose[df_noClose$lifeStatus %in% c(\"Return\"), ]\n",
    "\n",
    "df_day_return$lifeStatus <- factor(df_day_return$lifeStatus)\n",
    "\n",
    "\n",
    "df_day_return_histogram_combo <- df_day_return %>%\n",
    "  group_by(userIndex, day_of_week, lifeStatus) %>%\n",
    "  summarize(count_per_user = n())\n",
    "\n",
    "df_day_return_histogram_combo$day_of_week <- factor(df_day_return_histogram_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "df_day_break_return_histogram_combo <- rbind(df_day_break_histogram_combo, df_day_return_histogram_combo)\n",
    "\n",
    "options(repr.plot.width=18, repr.plot.height=8)\n",
    "\n",
    "# Create the histogram plot\n",
    "ggplot(df_day_break_return_histogram_combo, aes(x = day_of_week)) +\n",
    "  geom_boxplot( aes(y = count_per_user, fill=lifeStatus)) +\n",
    "  labs(x = \"Day of the Week\", y = \"Count per User\", fill=\"Status\") +\n",
    "  scale_color_manual(values = two_color) + \n",
    "  scale_y_continuous(expand=c(0,0), limits=c(0, 5))+\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), legend.position = \"bottom\", legend.box=\"horizontal\")\n",
    "\n",
    "\n",
    "# ggsave('figures/day_of_week_box.pdf', width=12, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "m = art(count_per_user ~ day_of_week*lifeStatus, data = df_day_break_return_histogram_combo)\n",
    "anova(m) # sig\n",
    "# art.con(m, ~day_of_week*lifeStatus) # sig Break-return and dropout-return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Apps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=10, repr.plot.height=7)\n",
    "\n",
    "df_app_categories[df_app_categories$app_mapped == 'News & Entertainment',]$app_mapped <- 'News &\\nEntertainment'\n",
    "\n",
    "df_app_categories$app_mapped <- factor(df_app_categories$app_mapped, levels=df_app_categories$app_mapped)\n",
    "\n",
    "# df_day_return_histogram_combo$day_of_week <- factor(df_day_return_histogram_combo$day_of_week, levels = days_ordered)\n",
    "\n",
    "# Plot the histogram\n",
    "app_category_plot <- ggplot(df_app_categories, aes(x = app_mapped, y = percentage)) +\n",
    "  geom_bar(stat = \"identity\", fill = two_color[1], color = 'black', alpha = 0.5) +\n",
    "  geom_text(aes(label = round(percentage,1)), vjust = -0.5, size=7.5) +\n",
    "  labs(x = \"Target App Category\", y = \"% of all App Open Attempts\") +\n",
    "#   scale_x_discrete(limits=c(NA,80), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,100), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=0, b=0, l=5, unit='pt'), axis.text.x = element_text(angle = 45, hjust = 1))\n",
    "\n",
    "app_category_plot\n",
    "\n",
    "# ggsave('figures/app_categories.pdf', width=10, height=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=10, repr.plot.height=7)\n",
    "\n",
    "df_app_usage_temp <- df_app_usage[0:10,]\n",
    "\n",
    "df_app_usage_temp$app <- c('Instagram', 'Twitter', 'YouTube', 'TikTok', 'Facebook', 'Reddit', 'Snapchat', 'Apollo (Reddit)', 'Telegram', 'LinkedIn')\n",
    "\n",
    "df_app_usage_temp$app <- factor(df_app_usage_temp$app, levels=df_app_usage_temp$app)\n",
    "\n",
    "\n",
    "\n",
    "# Plot the histogram\n",
    "app_plot <- ggplot(df_app_usage_temp, aes(x = app, y = percentage)) +\n",
    "  geom_bar(stat = \"identity\", fill = two_color[1], color = 'black', alpha = 0.5) +\n",
    "  geom_text(aes(label = round(percentage,1)), vjust = -0.5, size=7.5) +\n",
    "  labs(x = \"Target App\", y = \"% of all App Open Attempts\") +\n",
    "#   scale_x_discrete(limits=c(NA,80), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,100), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=0, b=0, l=5, unit='pt'), axis.text.x = element_text(angle = 45, hjust = 1))\n",
    "\n",
    "app_plot\n",
    "\n",
    "# ggsave('figures/apps_top10.pdf', width=10, height=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interventions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "intervention_user_counts <- df_noClose %>%\n",
    "  group_by(interventionType) %>%\n",
    "  summarize(user_count = n_distinct(userIndex)) %>%\n",
    "  ungroup()%>%\n",
    "  arrange(desc(user_count))\n",
    "\n",
    "intervention_user_counts <- intervention_user_counts[!is.na(intervention_user_counts$interventionType),]\n",
    "intervention_user_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "intervention_user_counts$interventionType <- c('Breathing Exercise', 'Minimal Breathing\\nExercise', 'Follow the Dot', 'Mirror', 'Black Screen', 'Rotate Phone')\n",
    "intervention_user_counts$interventionType <- factor(intervention_user_counts$interventionType, levels=intervention_user_counts$interventionType)\n",
    "\n",
    "options(repr.plot.width=20, repr.plot.height=7)\n",
    "\n",
    "# Plot the histogram\n",
    "interventions_plot <- ggplot(intervention_user_counts, aes(x = interventionType, y = user_count)) +\n",
    "  geom_bar(stat = \"identity\", fill = two_color[1], color = 'black', alpha = 0.5) +\n",
    "  geom_text(aes(label = user_count), vjust = -0.5, size=7.5) +\n",
    "  labs(x = \"Intervention Type\", y = \"Unique User Count\") +\n",
    "  # scale_x_discrete(limits=c(NA,80), expand=c(0,0),  breaks=pretty_breaks()) +\n",
    "  scale_y_continuous(limits=c(0,1010), expand=c(0,0), breaks=pretty_breaks()) +\n",
    "  theme_classic() +\n",
    "  theme(text = element_text(size=25), plot.margin = margin(t=10, r=0, b=0, l=0, unit='pt'), axis.text.x = element_text(angle = 45, hjust = 1))\n",
    "\n",
    "interventions_plot\n",
    "\n",
    "# ggsave('figures/interventions.pdf', width=20, height=7)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.0"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}