Twitter, Poisson Processes, and Conjugacy

Pulling Data from Twitter’s API

First I pulled data using the code below. For this code to work, your API keys will need to be specified in the setup above. To conceil my keys, I ran the commented code below earlier and saved the timeline object. The uncommented code loads my timeline and looks at the first element. The timeline is represented as a list of Status objects like the one output by the code below.

# # Twitter handel to pull data from
# handle = '@macstrelioff'
# 
# # get timeline
# timeline=api.GetUserTimeline(screen_name = handle, 
#                     count=200, # 200 is maximum 
#                     include_rts=True, 
#                     trim_user=True, 
#                     exclude_replies=False)
# 
# # save timeline object
# save_obj(timeline,'timeline_macstrelioff_20190406')

# load timeline object
timeline=load_obj('timeline_macstrelioff_20190406')
timeline[0] # most recent tweet status object

>>> Status(ID=1113860296458756097, ScreenName=None, Created=Thu Apr 04 17:46:03 +0000 2019, Text="@vboykis df.dropna(how='brute force') https://t.co/QOULUc5a0u")

Working with Twitter Status objects

I cover three ways to work with Status objects. 1. display the Status as a tweet! 2. Convert the Status to a dictionary and access values from keys 3. Access values directly as attributes of the Status

Embed a Status

First, many the Status attributes (created_at, favorite_count, text, …) can be cleanly displayed as a tweet embedded in a notebook. Below I create function that takes a username and tweet ID then, using Twitter’s embedding API, displayes the tweet as it would be seen on Twitter. (Note: this will only work properly if the Python kernel is trusted)

# for displaying tweets based on username and tweet_id
class disp_tweet(object):
    def __init__(self, user_name, tweet_id):
        # see: https://dev.twitter.com/web/embedded-tweets
        api = 'https://publish.twitter.com/oembed?url=https://twitter.com/'+ \
               user_name + '/status/' + tweet_id
        response  = requests.get(api)
        self.text = response.json()["html"]

    def _repr_html_(self):
        return self.text

disp_tweet(user_name='macstrelioff',tweet_id='981338927419109376')

>>> <__main__.disp_tweet object at 0x1a26d06710>

If run from a Jupyter notebook, this should embed a tweet as below;

My desk is covered in random papers. It is the support of a stationery distribution.

— mac strelioff (@macstrelioff) April 4, 2018

Convert a Status to a dict

Status objects have a bound method, .AsDict(), that will convert them to a Python dictionary. This way the structure of the information is easily seen. In the code below, I convert the first status to a dictionary and output it contents.

print(timeline[0].AsDict())

>>> {'created_at': 'Thu Apr 04 17:46:03 +0000 2019', 'favorite_count': 20, 'hashtags': [], 'id': 1113860296458756097, 'id_str': '1113860296458756097', 'in_reply_to_screen_name': 'vboykis', 'in_reply_to_status_id': 1113822568211996672, 'in_reply_to_user_id': 19304217, 'lang': 'da', 'media': [{'display_url': 'pic.twitter.com/QOULUc5a0u', 'expanded_url': 'https://twitter.com/macstrelioff/status/1113860296458756097/photo/1', 'id': 1113860285566046210, 'media_url': 'http://pbs.twimg.com/tweet_video_thumb/D3U6BzqUUAIwYEC.jpg', 'media_url_https': 'https://pbs.twimg.com/tweet_video_thumb/D3U6BzqUUAIwYEC.jpg', 'sizes': {'thumb': {'w': 150, 'h': 150, 'resize': 'crop'}, 'large': {'w': 250, 'h': 198, 'resize': 'fit'}, 'medium': {'w': 250, 'h': 198, 'resize': 'fit'}, 'small': {'w': 250, 'h': 198, 'resize': 'fit'}}, 'type': 'animated_gif', 'url': 'https://t.co/QOULUc5a0u', 'video_info': {'aspect_ratio': [125, 99], 'variants': [{'bitrate': 0, 'content_type': 'video/mp4', 'url': 'https://video.twimg.com/tweet_video/D3U6BzqUUAIwYEC.mp4'}]}}], 'retweet_count': 2, 'source': '<a href="http://twitter.com/download/android" rel="nofollow">Twitter for Android</a>', 'text': "@vboykis df.dropna(how='brute force') https://t.co/QOULUc5a0u", 'urls': [], 'user': {'id': 70255183, 'id_str': '70255183'}, 'user_mentions': [{'id': 19304217, 'id_str': '19304217', 'name': 'Vicki Boykis', 'screen_name': 'vboykis'}]}

{'created_at': 'Thu Apr 04 17:46:03 +0000 2019',
 'favorite_count': 20,
 'hashtags': [],
 'id': 1113860296458756097,
 'id_str': '1113860296458756097',
 'in_reply_to_screen_name': 'vboykis',
 'in_reply_to_status_id': 1113822568211996672,
 'in_reply_to_user_id': 19304217,
 'lang': 'da',
 'media': [{'display_url': 'pic.twitter.com/QOULUc5a0u',
   'expanded_url': 'https://twitter.com/macstrelioff/status/1113860296458756097/photo/1',
   'id': 1113860285566046210,
   'media_url': 'http://pbs.twimg.com/tweet_video_thumb/D3U6BzqUUAIwYEC.jpg',
   'media_url_https': 'https://pbs.twimg.com/tweet_video_thumb/D3U6BzqUUAIwYEC.jpg',
   'sizes': {'thumb': {'w': 150, 'h': 150, 'resize': 'crop'},
    'large': {'w': 250, 'h': 198, 'resize': 'fit'},
    'medium': {'w': 250, 'h': 198, 'resize': 'fit'},
    'small': {'w': 250, 'h': 198, 'resize': 'fit'}},
   'type': 'animated_gif',
   'url': 'https://t.co/QOULUc5a0u',
   'video_info': {'aspect_ratio': [125, 99],
    'variants': [{'bitrate': 0,
      'content_type': 'video/mp4',
      'url': 'https://video.twimg.com/tweet_video/D3U6BzqUUAIwYEC.mp4'}]}}],
 'retweet_count': 2,
 'source': '<a href="http://twitter.com/download/android" rel="nofollow">Twitter for Android</a>',
 'text': "@vboykis df.dropna(how='brute force') https://t.co/QOULUc5a0u",
 'urls': [],
 'user': {'id': 70255183, 'id_str': '70255183'},
 'user_mentions': [{'id': 19304217,
   'id_str': '19304217',
   'name': 'Vicki Boykis',
   'screen_name': 'vboykis'}]}

Access Status Attributes (How I accessed the data used below)

Since I’m interested in modeling expected number of tweets in a week, the most relevant attribute is the timestamps in the created_at attribute. These attributes can be accessed directly from the Status object. Below I make a list of the times at which each tweet was created and check the first element of that list;

# get list of time stamps
times = [pd.Timestamp(tweet.created_at) for tweet in timeline]
times.reverse() # sort s.t. times[0] is lowest, times[-1] is highest
times[0]

>>> Timestamp('2018-11-01 05:56:58+0000', tz='tzutc()')

This creates a list of timestamps. By default, times from the twitter API are localized to the UTC timezone. Below I convert these to my local time in California;

times = [time.tz_convert("America/Los_Angeles") for time in times]
times[0]

>>> Timestamp('2018-10-31 22:56:58-0700', tz='America/Los_Angeles')

Now we have a list of timestamps in local time! To get a sense of the duration over which this data spans, below I compute the time difference between the frist and last timestamp;

times[-1]-times[0]

>>> Timedelta('154 days 11:49:05')

Woah, almost 155 days of my twitter activity!

Assumptions

1. No more than one event can occur at a single point in time.
   • This can be violated when a user publishes a thread of multiple tweets at once. This can be fixed by recoding threads as a single status.
2. Independence: The interval lengths for each event are not influenced by any other event.
   • This can be violated if, instead of using Twitter’s thread option, a user ends a tweet with “…” to indicate that they will soon create another tweet. In this case, there are some tweets that imply a shorter interval before the next tweet.
3. Homogeneity: The distribution of intervals is the same throughout the entire process.
   • I probe this assumption in depth below, and it almost certainly violated.
   • There are methods for modeling inhomogeneous Poisson processes, but I ignore those here.

Specification and Properties

In this context tweet events are occurring across time. I index tweets with \(i\in\{1,...,N\}\), where \(N\) is the total number of tweets observed. Each tweet is created at a time, \(t_i\), and the next tweet is observed after an interval \(s_{i}\). That is, if tweet \(i\) is created at time \(t_{i-1}\) then tweet \(i+1\) is created at time \(t_{i}=t_{i-1}+s_{i}\). The interval between each tweet is \(s_i = t_{i}-t_{i-1} = (t_{i-1}+s_i)-t_{i-1}\). The assumptions of a Poisson process permit the following distributions for three interesting features of this scenario.

1. The distribution of time between events, \(s_i\), is exponential; \(s_i\sim Expo(\lambda) \Rightarrow p(s_i|\lambda) = \lambda e^{-\lambda s_i}\)
   • Here \(\lambda\) is a parameter that describes the tweet rate.
2. Given an interval of length \(s\), the distribution of the count of events in that interval, \(c|s\), is Poisson; \(c|s\sim Poisson(\lambda s) \Rightarrow p(c|s,\lambda) = \frac{(\lambda s)^{c}e^{-\lambda s}}{c!}\)
   • The count of events, \(c\), in a fixed interval depends both on the rate of the events, \(\lambda\), and the duration of the interval, \(s\).
3. The distribution of the total interval required for \(c\) events, \(s|c\), is gamma; \(s|c \sim Gamma(c,\lambda) \Rightarrow p(s|c,\lambda) = \frac{\lambda^c}{\Gamma(c)}(s)^{c-1}e^{-\lambda s}\)
   • The interval, \(s\), required for a fixed number of events depends on both the rate of the events, \(\lambda\), and the number of events, \(c\).

More information on these three kinds of distributions, and ways to implement them in Python, can be found in the scipy documentation on statistical functions. General information on each of these distributions can be found on the Wikipedia page for the exponential, Poisson, or gamma distribution. A key difference between the standard uses of these distributions and their roles in a Poisson process is that the rate parameter \(\lambda\) is also scaled by the duration of an interval \(s\) when constructing a distribution for the count of events in interval \(s\) (2, above) or the duration of the interval required for \(c\) events (3, above).

Checking the Homogeneity Assumption

The homogeneity assumption strictly requires that tweet rates are constant across time. This would generate data that are uniform across meaningful intervals such as time in a week or time in a day. To check homogeneity, below I convert the timestamps into the hour within a week, minute within a day, and minute within an hour, and plot tweet counts across these representations of time.

# convert to hours in a week
hour_of_week     = [t.weekday()*24+t.hour+t.minute/60 for t in times]
minute_of_day    = [t.hour*60+t.minute+t.second/60 for t in times]
minute_of_hour   = [t.minute+t.second/60 for t in times]
plt.figure(figsize=(10,4));
plt.hist(hour_of_week,bins=80);

plt.title("My Tweet Counts By Hour Of Week");
plt.ylabel("Count");
plt.xlabel("Hour In Week");

The histogram above indicates that there might be some hours of the week that have a higher rate than others. For example, I don’t seem to tweet much early on Sunday (hours 0-5), but I do seem to tweet a lot during the day on Sunday (around hours 6-20). Below I use rug plots, which represent a tweet event with a vertical line near the x-axis, and an imposed kernel density estimate, which is a continuous version of a histogram. I remake this plot in terms of hours within a week, minutes within a day, and minutes within an hour. If the tweet rate (\(\lambda\)) were homogeneous, then the kernel density estimate would be approximately flat.

def rug_plot_and_density(dat,bw,xlab,xlim):
    # rug plot + density
    plt.figure(figsize=(10,4))
    sns.distplot(dat, hist = False, kde = True, rug = True,
                 color = 'darkblue', 
                 kde_kws={'linewidth': 3,"bw":bw},
                 rug_kws={'color': 'black'})
    # formatting
    plt.title('Tweet Density By '+xlab)
    plt.xlabel(xlab)
    plt.ylabel('Kernel Density')
    plt.xlim(xlim);
    plt.show()

# data and plot formatting arguments
dats = (hour_of_week,minute_of_day,minute_of_hour)
xlabs=('Hour of Week','Minute of Day','Minute of Hour')
xlims=([0,24*7],[0,24*60],[0,60]);
bws  = (4,40,4)

# make plots
for dat,xlab,xlim,bw in zip(dats,xlabs,xlims,bws):
    rug_plot_and_density(dat=dat,bw=bw,xlab=xlab,xlim=xlim)

Are very sensitive to the choice of the bandwidth parameter that determines the window over which to aggregate events (similar to bin size when using a histogram). I’m using these plots to demonstrate possible violations of homogeneity that I would follow up on in a real analysis, but will not follow up on here.

From the top plot, there seems to be two patterns. First, a series of peaks and troughs that roughly correspond to daytime and night-time hours. I probably tweet with a higher frequency when I am awake, rather than asleep – meaning that \(\lambda\) may depend on time within a day. Second, a generally lower kernel density estimate during the middle than the edges. I may tweet more during the weekends (edge hours) than week days (middle hours).

From the middle plot that displays tweet frequencies by minutes within a day, there again seem to be two trends. First, I rarely tweet before minute 400 (around 6:40AM). Second, I have peaks around minute 600 (10:00AM), 1000 (4:40PM), and 1350 (10:00PM). This might be related to the times that I take a break from working. I generally take a break around 5:00PM, and usually take another break before bed around 9:00-10:00PM.

The bottom plot displays tweet frequencies by minutes within an hour. This seems more flat overall longer periods of time, but I may strangely tend to tweet more during the first half of hours.

Overall, there are many reasons that the homogeneity assumption may be violated. For cases like this, the tweet rate \(\lambda\) can be modeled as a function of time. However, to keep this example simple, I’ll ignore possible violations and proceede as if \(\lambda\) were a constant with respect to time.

Checking the exponential distribution of intervals

Let’s the distribution of the time between tweets, \(s_i\). If the assumptions of the Poisson process were satisfied, then the intervals between tweets would follow an exponential distribution. Below I compute the number of seconds between tweets and display each value as a black dash on the x-axis. I overlay a histogram, a kernel density, and a exponential density based on the observed mean interval.

# compute intervals between tweets
ss = [times[i]-times[i-1] for i in range(1,len(times))];
# convert from Timedelta to total time in seconds
ss = [si.total_seconds() for si in ss]

# histogram and density plot
plt.figure(figsize=(10,4))
sns.distplot(ss, hist = True, kde = True, rug = True,
             color = 'darkblue', bins=100,
             hist_kws={'color':[0,.7,.5,.5],'label':'Histogram'},
             kde_kws={'linewidth': 3,"bw":60*60,'label':'Kernel Density'},
             rug_kws={'color': 'black'})
plt.xlim([0,680000]);

plt.title('Distribution of intervals between tweets')
plt.xlabel('Seconds')
plt.ylabel('Frequency')

# overlay an exponential density
tmp_rate=np.mean(ss)
tmpx = np.linspace(0,680000,680000*5)
tmpy = expon.pdf(tmpx,scale=tmp_rate)
plt.plot(tmpx,tmpy,color=[.7,0,0,1],linewidth=3,label='Exponential Density');
# add legend
plt.legend();

Based on the relative heights of the kernel density and the exponential density, there seem to be more short intervals, fewer moderate length intervals, and more long intervals relative to the exponential distribution. The mean and standard deviation of an exponential distribution should be the same value. Below I check the standard deviation and mean of the observed intervals.

# check standard deviation
np.std(ss),np.mean(ss),np.std(ss)/np.mean(ss)

>>> (90150.9122414953, 67076.1055276382, 1.3440093388300445)

The observed standard deviation is about 1.34 times larger (variance is about 1.80 times larger) than it would be if the data were exponentially distributed with the observed mean. While this could be accounted for with an overdispersion parameter, I will ignore this issue here for the sake of having a simple and fast online model.

Inference On Tweet Rate \(\lambda\) Over Time

In the code below, I compute the maximum a posteriori (MAP) estimate, or posterior mode, across time. I plot this across time with the 97.5\(^{th}\) and 2.5\(^{th}\) percentiles as a shaded region representing the 95% credible interval.

# get MAP lambda
lambda_map   = [(alpha-1)/beta for alpha,beta in zip(alphas,betas)]
# get CI upper and lower bounds
lambda_upper = [gamma.ppf(.975,a=alpha,scale=1/beta) for alpha,beta in zip(alphas,betas)]
lambda_lower = [gamma.ppf(.025,a=alpha,scale=1/beta) for alpha,beta in zip(alphas,betas)]
# plot over time
plt.figure(figsize=(10,4))
plt.plot(times[1:],lambda_map,color=[0,0,1],label='MAP');
plt.fill_between(times[1:], lambda_lower, lambda_upper,color=[0,0,1,.1],label='95% CI')
#sns.distplot(times[1:], hist = False, kde = False, rug = True,
#             color = 'darkblue', 
#             rug_kws={'color': 'black'});
plt.ylabel('Posterior on Lambda (tweets/day)')
plt.xlabel('Time')
plt.title('Summary of the posterior over $\lambda$ across time')
plt.legend();

The posterior seems to tighten dramatically around the mode during the first month, and the mode seems to stabalize after about that much time. However, this method seems slow to adjust to the slower rate of tweets in 2019. One way to address this might be to add weights to the parameter updates such that \(\alpha\) and \(\beta\) are more sensitive to recent data than past data. These weights (or stepsize or learning rates) can be based on the surprisal, or likelihood, of a new tweet under the posterior. If the new interval was well anticipated, then little updating is needed, but if the new interval was very surprising or unlikely, then a sharpe change in \(\alpha\) and \(\beta\) might be warrented.