Fonterra Social Media Analysis

Overview

This report analyses 9,001 Fonterra social media posts from six platforms (Facebook, Instagram, LinkedIn, TikTok, Twitter/X and YouTube), published between 2010 and 2026.

The analysis has three parts:

• Data preparation: cleaning the timestamps, media labels and missing engagement counts, and preparing the post text for analysis.
• Text analysis: finding the most frequent words and using topic modelling (LDA, 6 topics) to see what Fonterra posts about.
• Engagement analysis: comparing likes, comments and shares across platforms, media types and posting patterns.

Install and load packages and import and preview data

# Set consistent options for tables and figures
knitr::opts_chunk$set(fig.width = 8, fig.height = 5, echo = TRUE)

# Install and load packages
if(!require("pacman"))install.packages("pacman")

## Loading required package: pacman

pacman::p_load(readr, dplyr, tidyr, stringr, lubridate, purrr, ggplot2, gridExtra, GGally, corrplot,forcats,rstatix,DescTools, skimr, finalfit)
pacman::p_load(tidytext, textstem, wordcloud, textdata, topicmodels, reshape2, 
               knitr)

# Import the Fonterra social media data
raw_post_data <- read_csv(
  "fonterra.csv",
  na = c("", "NA"),
  show_col_types = FALSE
)

# Preview the structure and first rows of the data
glimpse(raw_post_data)

## Rows: 9,001
## Columns: 9
## $ platform       <chr> "Twitter/X", "Twitter/X", "Twitter/X", "Twitter/X", "Tw…
## $ url            <chr> "https://x.com/Fonterra/status/1891573780557988062", "h…
## $ createdAt      <chr> "2025/02/17 19:41:47", "2024/12/05 22:23:33", "2024/12/…
## $ text           <chr> "Today we announced new incentives to help farms reduce…
## $ media          <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ commentCount   <dbl> 66, 1, 1, 303, 8, 2, 2, 0, 2, 0, 1, 1, 1, 0, 0, 1, 0, 0…
## $ shareCount     <dbl> 2, 0, 0, 102, 3, 0, 0, 0, 2, 0, 2, 0, 0, 0, 0, 1, 0, 0,…
## $ likeCount      <dbl> 15, 0, 0, 704, 16, 1, 6, 0, 19, 0, 16, 1, 3, 18, 2, 14,…
## $ videoplayCount <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…

head(raw_post_data)

## # A tibble: 6 × 9
##   platform  url          createdAt text  media commentCount shareCount likeCount
##   <chr>     <chr>        <chr>     <chr> <chr>        <dbl>      <dbl>     <dbl>
## 1 Twitter/X https://x.c… 2025/02/… "Tod… <NA>            66          2        15
## 2 Twitter/X https://x.c… 2024/12/… "@ah… <NA>             1          0         0
## 3 Twitter/X https://x.c… 2024/12/… "@ah… <NA>             1          0         0
## 4 Twitter/X https://x.c… 2024/12/… "At … <NA>           303        102       704
## 5 Twitter/X https://x.c… 2024/11/… "We’… <NA>             8          3        16
## 6 Twitter/X https://x.c… 2024/10/… "Tod… <NA>             2          0         1
## # ℹ 1 more variable: videoplayCount <dbl>

Data preparation

Non-text data

There are some criteria taken into account during the non-text data cleaning:

• createdAt variable was changed to UTC date-time, since according to the instruction, all the timestamps were collected in UTC time zone.
• There are some timestamps with 00:00:00, which can be interpreted as date only without time information. Therefore, time will not be considered, only date-related features (year, month, weekday).
• Likes and “digs” on TikTok are regarded as similar. Comments (except Twitter/X) and replies (Twitter/X) are also considered to be similar.
• media feature contains different capitalisation (video, Video), so it should be standardised (Image, Video, Sidecar). Empty value was assigned Not recorded, rather than “no media”, since media type cannot be determined.
• High levels of engagement are maintained, since they probably are real examples of viral posts rather than outliers; rather than deleting them, medians and non-parametric tests are employed later on.
• Posts without text are maintained during engagement analysis but deleted only when text analysis begins, since there will be nothing to tokenize then.
• The missing engagements remain missing as opposed to using 0s since certain sites such as Instagram and YouTube do not provide the sharing metrics.
• There are no duplicates in the form of identical rows or URLs. However, there is one instance where two posts have the same date and text for Facebook and Instagram (January 11, 2019, 06:00:24). This is not an error but cross-posting.

# confirm all expected columns are present before cleaning
necessary_cols <- c("platform", "url", "createdAt", "text", "media",
                  "commentCount", "shareCount", "likeCount", "videoplayCount")
stopifnot(all(necessary_cols %in% names(raw_post_data)))

# trim whitespace, parse the timestamp, and standardise the media labels
posts <- raw_post_data |>
  mutate(
    platform   = str_squish(platform),
    url        = str_squish(url),
    createdAt  = ymd_hms(createdAt, tz = "UTC", quiet = TRUE),
    post_media = case_when(
      is.na(media)                             ~ "Not recorded",
      str_to_lower(media) %in% c("photo", "image") ~ "Image",
      str_to_lower(media) == "video"            ~ "Video",
      str_to_lower(media) == "sidecar"          ~ "Sidecar",
      TRUE ~ str_to_title(media)
    )
  )

# confirm every timestamp parsed successfully
stopifnot(sum(is.na(posts$createdAt)) == 0)

# check for exact duplicate rows and duplicate URLs
tibble(
  check  = c("Exact duplicate rows", "Duplicate URLs"),
  number = c(sum(duplicated(posts)), sum(duplicated(posts$url)))
) |>
  kable(caption = "Duplicate checks")

Duplicate checks

check	number
Exact duplicate rows	0
Duplicate URLs	0

# check for identical text posted at the exact same time, and which platforms
posts |>
  filter(!is.na(text)) |>
  group_by(createdAt, text) |>
  filter(n() > 1) |>
  summarise(platforms = paste(platform, collapse = " & "), records = n(), .groups = "drop") |>
  select(createdAt, platforms, records) |>
  kable(caption = "Same-time, same-text records (kept as likely cross-posting)")

Same-time, same-text records (kept as likely cross-posting)

createdAt	platforms	records
2019-01-11 06:00:24	Facebook & Instagram	2

# missing values by column
posts |>
  summarise(across(everything(), ~ sum(is.na(.x)))) |>
  pivot_longer(everything(), names_to = "column", values_to = "missing") |>
  mutate(pct_missing = round(100 * missing / nrow(posts), 1)) |>
  kable(caption = "Missing values by column")

Missing values by column

column	missing	pct_missing
platform	0	0.0
url	0	0.0
createdAt	0	0.0
text	82	0.9
media	6994	77.7
commentCount	293	3.3
shareCount	951	10.6
likeCount	0	0.0
videoplayCount	8115	90.2
post_media	0	0.0

# check whether missing engagement counts are platform-specific
posts |>
  group_by(platform) |>
  summarise(
    missing_comments = sum(is.na(commentCount)),
    missing_shares    = sum(is.na(shareCount)),
    missing_plays     = sum(is.na(videoplayCount)),
    .groups = "drop"
  ) |>
  kable(caption = "Missing engagement counts by platform")

Missing engagement counts by platform

platform	missing_comments	missing_shares	missing_plays
Facebook	269	0	1115
Instagram	0	519	454
LinkedIn	24	37	326
TikTok	0	0	0
Twitter/X	0	0	6220
YouTube	0	395	0

# check for impossible negative counts
posts |>
  summarise(across(c(commentCount, shareCount, likeCount, videoplayCount),
                    ~ sum(.x < 0, na.rm = TRUE))) |>
  pivot_longer(everything(), names_to = "column", values_to = "negative_values") |>
  kable(caption = "Any impossible negative counts?")

Any impossible negative counts?

column	negative_values
commentCount	0
shareCount	0
likeCount	0
videoplayCount	0

New variables created

Why each new variable is useful:

• post_id gives each post a unique identifier which is required for the text analysis and topic model.
• post_date, post_year, post_month, post_yearmonth and post_weekday break the timestamp into components, so posting patterns can be examined without altering the original createdAt.
• post_media is a tidy content-format label that can be compared against engagement.
• word_count and char_count measure caption length, so length can be tested against engagement.
• total_reaction sums whichever of likes/comments/shares a post has data for. It is a useful descriptive number but not the main cross-platform outcome, since shares are structurally missing on some platforms.
• has_play_data flags the subset of posts where a plays-vs-likes comparison is possible.
• has_hashtag flags whether a post used a # - a two-level variable used for a t-test/Mann-Whitney comparison later.

#new_variables
posts <- posts |>
  mutate(
    post_id        = row_number(),
    post_date      = as.Date(createdAt),
    post_year      = year(createdAt),
    post_month     = month(createdAt, label = TRUE, abbr = TRUE),
    post_yearmonth = floor_date(createdAt, unit = "month"),
    post_weekday   = wday(createdAt, label = TRUE, abbr = FALSE, week_start = 1),
    word_count     = str_count(coalesce(text, ""), boundary("word")),
    char_count     = str_length(coalesce(text, "")),
    total_reaction = likeCount + coalesce(commentCount, 0) + coalesce(shareCount, 0),
    has_play_data  = !is.na(videoplayCount),
    has_hashtag    = str_detect(coalesce(text, ""), "#\\w+")
  )

# preview the new variables
posts |>
  select(post_id, platform, createdAt, post_year, post_month, post_weekday,
         post_media, word_count, total_reaction, has_play_data) |>
  head(10) |>
  kable(caption = "New variables - first 10 rows")

New variables - first 10 rows

post_id	platform	createdAt	post_year	post_month	post_weekday	post_media	word_count	total_reaction	has_play_data
1	Twitter/X	2025-02-17 19:41:47	2025	Feb	Monday	Not recorded	54	83	FALSE
2	Twitter/X	2024-12-05 22:23:33	2024	Dec	Thursday	Not recorded	23	1	FALSE
3	Twitter/X	2024-12-05 21:05:34	2024	Dec	Thursday	Not recorded	27	1	FALSE
4	Twitter/X	2024-12-01 23:04:24	2024	Dec	Sunday	Not recorded	24	1109	FALSE
5	Twitter/X	2024-11-10 19:39:51	2024	Nov	Sunday	Not recorded	44	27	FALSE
6	Twitter/X	2024-10-10 19:36:23	2024	Oct	Thursday	Not recorded	46	3	FALSE
7	Twitter/X	2024-10-09 19:38:11	2024	Oct	Wednesday	Not recorded	51	8	FALSE
8	Twitter/X	2024-10-09 19:22:00	2024	Oct	Wednesday	Not recorded	54	0	FALSE
9	Twitter/X	2024-09-29 19:43:38	2024	Sep	Sunday	Not recorded	39	23	FALSE
10	Twitter/X	2024-09-28 04:00:18	2024	Sep	Saturday	Not recorded	63	0	FALSE

The posts run from 08 November 2010 to 28 July 2026. Because likes/comments/shares are cumulative counts rather than a rate over a fixed window, older posts have had more time to accumulate engagement than newer ones.

Text data

Posts with no text content are excluded, as there is no data to tokenize. Links, @mentions, HTML artifacts, and the # symbol are stripped before tokenization takes place, although the text within the # symbol is usually kept, as it relates to the topic. Stop words, along with a few special dataset filler words, are excluded.Curly apostrophes have been converted into straight apostrophes to ensure that contractions like we’re are excluded from the list of stop words. No lemmatisation was performed; therefore, different forms of a word are considered independently, for example, award and awards.

# custom stop words: generic filler that adds no subject-matter meaning
additional_stopwords <- tibble(
  word = c("fonterra", "amp", "https", "http", "www", "com", "co",
           "nz", "new", "zealand", "uh", "um", "yeah", "really",
           "just", "know", "like", "good", "great", "thanks", "thank",
           "hi", "hey", "cheers", "today", "year", "years", "day",
           "time", "going", "make", "right", "week", "look", "think",
           "got", "need", "want", "lot", "share", "come", "way",
           "doing", "sure", "things", "bit", "read", "link")
)

# remove links, @mentions and HTML artefacts, lower-case, then tidy whitespace
posts_with_text <- posts |>
  filter(!is.na(text), str_squish(text) != "") |>
  mutate(
    text_clean = text |>
      str_to_lower() |>
      str_replace_all("[’‘]", "'") |>
      str_replace_all("https?://\\S+|www\\.\\S+", " ") |>
      str_replace_all("@\\w+", " ") |>
      str_replace_all("&amp;", " and ") |>
      str_replace_all("#", " ") |>
      str_squish()
  ) |>
  select(post_id, platform, createdAt, text_clean)

# tokenise into one word per row (unnest_tokens() lower-cases and strips
# punctuation automatically)
post_words_raw <- posts_with_text |>
  unnest_tokens(output = word, input = text_clean, token = "words")

# remove standard stop words, the custom list above, and any token with a digit
post_words <- post_words_raw |>
  filter(!word %in% stop_words$word,
         !word %in% additional_stopwords$word,
         !str_detect(word, "[0-9]"),
         str_length(word) >= 3)

# how many posts were available for text analysis
tibble(
  measure = c("All posts", "Posts with usable text", "Posts dropped (no text)"),
  value = c(nrow(posts), n_distinct(posts_with_text$post_id),
            nrow(posts) - n_distinct(posts_with_text$post_id))
) |>
  kable(caption = "Text preparation summary")

Text preparation summary

measure	value
All posts	9001
Posts with usable text	8919
Posts dropped (no text)	82

Data-informed insights and recommendations

Topics content creators have been posting about

Term frequency shows which individual words appear most often; the LDA topic model groups words that tend to occur together into broader context. The two are complementary.

# count each word, then keep the 20 most frequent
word_counts <- post_words |> count(word, sort = TRUE, name = "frequency")
top_words <- word_counts |> slice_max(n = 20, order_by = frequency)

kable(top_words, caption = "Top 20 words in Fonterra's social media content")

Top 20 words in Fonterra’s social media content

word	frequency
milk	2127
dairy	1512
farmers	1492
farm	1213
team	1147
people	767
food	661
world	604
business	575
support	537
products	518
water	502
cheese	490
site	485
price	450
learn	411
global	404
tanker	378
china	376
quality	370

# bar chart of the top 20 words, shaded by frequency
top_words |>
  mutate(word = fct_reorder(word, frequency)) |>
  ggplot(aes(x = word, y = frequency, fill = frequency)) +
  geom_col() +
  coord_flip() +
  scale_fill_gradient(low = "#d9f2d0", high = "#1b4332", guide = "none") +
  labs(title = "Top 20 words across Fonterra's social media content",
       subtitle = "How often each word shows up across every platform",
       x = NULL, y = "Number of times used") +
  theme_minimal(base_size = 12)

The most common words in this word list include milk, dairy, farmers, and farm; then come team, people, food, world, business, support, products, and water. Some other common words in this list are cheese, price, China, global, and tanker. The contents of this list are not confined only to marketing of products but include many other topics.

# drop words that are too rare (fewer than 5 posts) or too common (over half
# of all posts) - both extremes make topics harder to interpret
word_doc_counts <- post_words |>
  distinct(post_id, word) |>
  count(word, name = "n_docs")

n_text_posts <- n_distinct(post_words$post_id)

common_closer <- word_doc_counts |>
  filter(n_docs >= 5, n_docs <= 0.50 * n_text_posts)

input_topic_words <- post_words |>
  inner_join(common_closer, by = "word")

# build the document-term matrix that LDA() requires
topic_dtm <- input_topic_words |>
  count(post_id, word, name = "n") |>
  cast_dtm(document = post_id, term = word, value = n)

# fit a six-topic model (k = 6, as suggested in the brief); the seed in
# control= makes it reproducible
topic_model <- LDA(topic_dtm, k = 6, control = list(seed = 1234))

# quick list of the top words per topic (no probabilities attached)
get_terms(topic_model, 10)

##       Topic 1    Topic 2     Topic 3       Topic 4           Topic 5      
##  [1,] "milk"     "milk"      "support"     "dairy"           "farm"       
##  [2,] "farmers"  "dairy"     "school"      "tanker"          "farmers"    
##  [3,] "price"    "world"     "community"   "awards"          "water"      
##  [4,] "business" "food"      "kids"        "farmers"         "cows"       
##  [5,] "john"     "products"  "farmers"     "award"           "music"      
##  [6,] "million"  "anchor"    "farm"        "industry"        "farming"    
##  [7,] "china"    "cream"     "milk"        "congratulations" "dairy"      
##  [8,] "strong"   "nutrition" "team"        "agchatnz"        "environment"
##  [9,] "forecast" "cheese"    "people"      "team"            "people"     
## [10,] "food"     "protein"   "communities" "māori"           "quality"    
##       Topic 6    
##  [1,] "team"     
##  [2,] "site"     
##  [3,] "cheese"   
##  [4,] "hope"     
##  [5,] "check"    
##  [6,] "event"    
##  [7,] "christmas"
##  [8,] "gdt"      
##  [9,] "listen"   
## [10,] "global"

# word-topic probabilities, needed for the chart below
topic_words <- tidy(topic_model, matrix = "beta") |>
  group_by(topic) |>
  slice_max(order_by = beta, n = 10) |>
  ungroup()

topic_words |>
  arrange(topic, desc(beta)) |>
  group_by(topic) |>
  summarise(top_words = paste(term, collapse = ", "), .groups = "drop") |>
  kable(caption = "Top words for each of the six topics")

Top words for each of the six topics

topic	top_words
1	milk, farmers, price, business, john, million, china, strong, forecast, food
2	milk, dairy, world, food, products, anchor, cream, nutrition, cheese, protein
3	support, school, community, kids, farmers, farm, milk, team, people, communities
4	dairy, tanker, awards, farmers, award, industry, congratulations, agchatnz, team, māori
5	farm, farmers, water, cows, music, farming, dairy, environment, people, quality
6	team, site, cheese, hope, check, event, christmas, gdt, listen, global

# bar chart of the top words per topic, one colour per topic
topic_words |>
  ggplot(aes(x = beta, y = reorder_within(x = term, by = beta, within = topic),
             fill = factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free") +
  scale_y_reordered() +
  scale_fill_brewer(palette = "Set2") +
  labs(subtitle = "Top 10 highest-probability words per LDA topic",
       x = "Word-topic probability (beta)", y = NULL) +
  theme_minimal(base_size = 10)

# ---- Word cloud of the 80 most frequent words ----
set.seed(2026)
wordcloud(words = word_counts$word, freq = word_counts$frequency, max.words = 80, random.order = FALSE,
          scale = c(3.5, 0.6), colors = RColorBrewer::brewer.pal(8, "Dark2"))

# assign each post its single most likely topic, with the model's confidence
post_topics <- tidy(topic_model, matrix = "gamma") |>
  group_by(document) |>
  slice_max(gamma, n = 1, with_ties = FALSE) |>
  ungroup() |>
  transmute(post_id = as.integer(as.character(document)), topic, topic_confidence = gamma)

# share of posts falling into each topic
post_topics |>
  count(topic, name = "posts") |>
  mutate(pct = round(100 * posts / sum(posts), 1)) |>
  kable(caption = "Share of posts falling into each topic")

Share of posts falling into each topic

topic	posts	pct
1	1087	12.4
2	1471	16.8
3	1484	16.9
4	1535	17.5
5	1387	15.8
6	1803	20.6

Commentary on topic insights

The following six themes have been developed based on the most frequent words of each topic (numbers are random, hence description needs to be confirmed using the table “Top words for each of the six topics” on each knit). The first topic deals with business and pricing (price, million, China, strong, forecast); the second topic deals with products and nutrition (products, cream, nutrition, cheese, protein); the third topic deals with community and schools (support, school, community, kids, communities); the fourth topic deals with awards and industry recognition (awards, congratulations, industry, agchatnz) and captures milk tanker posts; the fifth topic deals with farming and environment (farm, water, cows, farming, environment); the sixth topic deals with team, events and site activity (team, site, event, Christmas, GDT, global) and is the biggest topic capturing 20.6% of the posts. The sixth topic is rather indistinctive, and words like “music” (topic 5) and “john” (topic 1) are ambiguous. Topic sizes are fairly balanced (12.4% to 20.6%). In this dataset, the most prominent topics extend beyond direct product-related content and include business, community, industry recognition, farming and events.

Level of user engagement and factors that influence user engagement

The likeCount can be applied to every single post, and therefore, it is considered a main dependent variable for cross-platform comparison. The comments, shares, and videos played are counted where possible; however, they cannot be considered a directly comparable variable due to lack of data on some platforms.

# Summarise engagement statistics by platform
safe_median <- function(x) if (all(is.na(x))) NA_real_ else median(x, na.rm = TRUE)

platform_engagement <- posts |>
  group_by(platform) |>
  summarise(
    posts           = n(),
    median_likes    = median(likeCount, na.rm = TRUE),
    iqr_likes       = IQR(likeCount, na.rm = TRUE),
    median_comments = safe_median(commentCount),
    median_shares   = safe_median(shareCount),
    median_plays    = safe_median(videoplayCount),
    .groups = "drop"
  )

platform_engagement |>
  kable(digits = 1, caption = "Engagement by platform - descriptive statistics")

Engagement by platform - descriptive statistics

platform	posts	median_likes	iqr_likes	median_comments	median_shares	median_plays
Facebook	1433	61	134.0	7	4	10172.5
Instagram	519	132	117.5	2	NA	3690.0
LinkedIn	326	141	141.2	5	4	NA
TikTok	108	7	8.2	0	2	2969.5
Twitter/X	6220	1	5.0	0	0	NA
YouTube	395	0	4.0	0	NA	505.0

# distribution of likes per post (log scale) - shows the heavy right skew
posts |>
  ggplot(aes(x = likeCount + 1)) +
  geom_histogram(bins = 40, fill = "#225ea8") +
  scale_x_log10() +
  labs(title = "Distribution of likes per post",
       subtitle = "Log scale; most posts receive few likes",
       x = "Likes + 1 (log scale)", y = "Number of posts") +
  theme_minimal(base_size = 11)

Most posts receive very few likes, while a small number receive many, so the distribution is strongly right-skewed. This is why medians and non-parametric tests are used.

# posting volume over time by platform
posts |>
  count(post_year, platform) |>
  ggplot(aes(x = post_year, y = n, fill = platform)) +
  geom_col() +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Number of posts per year by platform",
       x = "Year", y = "Posts", fill = "Platform") +
  theme_minimal(base_size = 10)

Twitter/X accounts for most of the posts in the dataset (6,220 of 9,001), so overall posting volume is largely driven by this platform.

#plot_engagement_by_platform
# boxplot of engagement by platform; log10(1 + likes) keeps a handful of
# viral posts from squashing the rest of the distribution
posts |>
  ggplot(aes(x = fct_reorder(platform, likeCount, .fun = median), y = log10(likeCount + 1),
             fill = platform)) +
  geom_boxplot(outlier.alpha = 0.25, show.legend = FALSE) +
  coord_flip() +
  scale_fill_brewer(palette = "Blues") +
  labs(title = "Engagement by platform",
       subtitle = "log10(1 + likes) per post, so a few viral posts don't hide everything else",
       x = NULL, y = "log10(1 + likes)") +
  theme_minimal(base_size = 11)

# boxplot of engagement by media type - only Facebook, Instagram and
# LinkedIn vary their media type in this dataset
posts |>
  filter(platform %in% c("Facebook", "Instagram", "LinkedIn")) |>
  ggplot(aes(x = post_media, y = log10(likeCount + 1), fill = post_media)) +
  geom_boxplot(outlier.alpha = 0.25) +
  facet_wrap(~ platform) +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = "Engagement by media type",
       subtitle = "Not recorded vs. image vs. video vs. sidecar, log10(1 + likes)",
       x = NULL, y = "log10(1 + likes)", fill = "Media type") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

#plot_length_vs_engagement
# scatter plot of caption length vs. likes; both axes logged since both
# variables are right-skewed
posts |>
  ggplot(aes(x = word_count + 1, y = likeCount + 1, colour = platform)) +
  geom_point(alpha = 0.35) +
  scale_x_log10() +
  scale_y_log10() +
  scale_colour_brewer(palette = "Set1") +
  labs(title = "Post length vs. engagement",
       subtitle = "Caption word count against likes (log-log scale), by platform",
       x = "Words in the post (log scale)", y = "Likes (log scale)", colour = "Platform") +
  theme_minimal(base_size = 11)

# boxplot of engagement by topic, linking Part 3(a)'s topics to Part 3(b)'s
# engagement question
posts |>
  left_join(post_topics, by = "post_id") |>
  filter(!is.na(topic)) |>
  ggplot(aes(x = factor(topic), y = log10(likeCount + 1), fill = factor(topic))) +
  geom_boxplot(outlier.alpha = 0.25, show.legend = FALSE) +
  scale_fill_brewer(palette = "Blues") +
  labs(title = "Engagement by topic",
       subtitle = "Does what a post is about relate to how well it does? (log10(1 + likes))",
       x = "Topic", y = "log10(1 + likes)") +
  theme_minimal(base_size = 11)

Prior to selecting any tests, the assumptions of normality (using Shapiro-Wilk test on each platform, but reducing sample size to 5,000 for those platforms which have more data since shapiro.test() can handle maximum 5,000 observations) and equality of variance (using rstatix::levene_test()) have been verified for the platform test while the assumption of Levene’s test has been verified for the hashtag test.

#bivariate_setup
# posts joined to their topic, and posts with video-play data - both reused below
posts_with_topic  <- posts |> left_join(post_topics, by = "post_id") |> filter(!is.na(topic)) |>
  mutate(topic = factor(topic))
videos_with_plays <- posts |> filter(!is.na(videoplayCount))

# log(1 + likes) column used for the normality/variance checks and the ANOVA
posts_log <- posts |> mutate(log_likes = log10(likeCount + 1))

# Shapiro-Wilk normality check, per platform (sampled to 5,000 rows where
# needed, since Twitter/X alone has over 6,000 posts)
posts_log |>
  group_by(platform) |>
  group_map(~ {
    vals <- .x$log_likes
    if (length(vals) > 5000) vals <- sample(vals, 5000)
    test <- shapiro.test(vals)
    tibble(platform = unique(.y$platform), n = length(.x$log_likes),
           W = unname(test$statistic), p_value = test$p.value)
  }) |>
  bind_rows() |>
  kable(digits = 3, caption = "Shapiro-Wilk normality check, log(1 + likes) by platform")

Shapiro-Wilk normality check, log(1 + likes) by platform

platform	n	W	p_value
Facebook	1433	0.998	0.066
Instagram	519	0.974	0.000
LinkedIn	326	0.992	0.065
TikTok	108	0.880	0.000
Twitter/X	6220	0.868	0.000
YouTube	395	0.753	0.000

# homogeneity-of-variance check
posts_log |> mutate(platform = factor(platform)) |> levene_test(log_likes ~ platform)

## # A tibble: 1 × 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1     5  8995      71.5 1.25e-73

Facebook and LinkedIn approach normality (p > .05), whereas Instagram, TikTok, Twitter/X, and YouTube deviate from normal (p < .001), with significant differences between variances for the different platforms (Levene’s test, p < .001). Hence, the Kruskal-Wallis outcome is considered the primary finding, with ANOVA included as well.

# one-way ANOVA (parametric) with Tukey's HSD post-hoc
anova_platform <- aov(log_likes ~ platform, data = posts_log)
summary(anova_platform)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## platform       5   3754   750.9    3257 <2e-16 ***
## Residuals   8995   2074     0.2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

as.data.frame(TukeyHSD(anova_platform)$platform) |>
  round(3) |>
  kable(caption = "Tukey's HSD - pairwise platform differences")

Tukey’s HSD - pairwise platform differences

	diff	lwr	upr	p adj
Instagram-Facebook	0.278	0.208	0.348	0.000
LinkedIn-Facebook	0.346	0.262	0.430	0.000
TikTok-Facebook	-0.862	-0.999	-0.725	0.000
Twitter/X-Facebook	-1.362	-1.403	-1.322	0.000
YouTube-Facebook	-1.409	-1.486	-1.331	0.000
LinkedIn-Instagram	0.068	-0.029	0.164	0.344
TikTok-Instagram	-1.140	-1.285	-0.995	0.000
Twitter/X-Instagram	-1.641	-1.703	-1.578	0.000
YouTube-Instagram	-1.687	-1.778	-1.595	0.000
TikTok-LinkedIn	-1.208	-1.360	-1.056	0.000
Twitter/X-LinkedIn	-1.708	-1.786	-1.631	0.000
YouTube-LinkedIn	-1.755	-1.857	-1.652	0.000
Twitter/X-TikTok	-0.500	-0.633	-0.368	0.000
YouTube-TikTok	-0.547	-0.695	-0.398	0.000
YouTube-Twitter/X	-0.046	-0.117	0.025	0.430

# Kruskal-Wallis (non-parametric equivalent) with Dunn's test as post-hoc
platform_test <- kruskal.test(likeCount ~ platform, data = posts)
platform_test

## 
##  Kruskal-Wallis rank sum test
## 
## data:  likeCount by platform
## Kruskal-Wallis chi-squared = 4728.4, df = 5, p-value < 2.2e-16

dunn_test(posts, likeCount ~ platform, p.adjust.method = "bonferroni") |>
  select(group1, group2, statistic, p.adj, p.adj.signif)

## # A tibble: 15 × 5
##    group1    group2    statistic     p.adj p.adj.signif
##    <chr>     <chr>         <dbl>     <dbl> <chr>       
##  1 Facebook  Instagram     4.85  1.84e-  5 ****        
##  2 Facebook  LinkedIn      4.55  7.92e-  5 ****        
##  3 Facebook  TikTok       -7.87  5.44e- 14 ****        
##  4 Facebook  Twitter/X   -53.1   0         ****        
##  5 Facebook  YouTube     -29.9   3.64e-195 ****        
##  6 Instagram LinkedIn      0.437 1   e+  0 ns          
##  7 Instagram TikTok       -9.77  2.22e- 21 ****        
##  8 Instagram Twitter/X   -39.5   0         ****        
##  9 Instagram YouTube     -29.2   8.40e-186 ****        
## 10 LinkedIn  TikTok       -9.59  1.36e- 20 ****        
## 11 LinkedIn  Twitter/X   -32.3   1.53e-227 ****        
## 12 LinkedIn  YouTube     -26.4   7.79e-153 ****        
## 13 TikTok    Twitter/X    -7.93  3.17e- 14 ****        
## 14 TikTok    YouTube      -8.41  5.90e- 16 ****        
## 15 Twitter/X YouTube      -2.77  8.46e-  2 ns

# chi-squared test: is media type linked to being above/below the platform's
# own median engagement? Restricted to Facebook, Instagram and LinkedIn,
# the only platforms where media type varies
posts_media <- posts |>
  filter(platform %in% c("Facebook", "Instagram", "LinkedIn")) |>
  group_by(platform) |>
  mutate(high_engagement = likeCount > median(likeCount)) |>
  ungroup()

chi_table1 <- table(posts_media$post_media, posts_media$high_engagement)
chi_table1

##               
##                FALSE TRUE
##   Image          434  496
##   Not recorded   456  318
##   Sidecar         46   49
##   Video          210  269

chisq.test(chi_table1)$expected   # check: all expected counts should be at least 5

##               
##                    FALSE      TRUE
##   Image        467.85777 462.14223
##   Not recorded 389.37840 384.62160
##   Sidecar       47.79192  47.20808
##   Video        240.97191 238.02809

chisq.test(chi_table1)

## 
##  Pearson's Chi-squared test
## 
## data:  chi_table1
## X-squared = 36.015, df = 3, p-value = 7.433e-08

chisq.test(chi_table1)$stdres

##               
##                     FALSE       TRUE
##   Image        -2.8866025  2.8866025
##   Not recorded  5.8943471 -5.8943471
##   Sidecar      -0.3756174  0.3756174
##   Video        -3.1849282  3.1849282

# hashtag use and median likes, by platform
posts |>
  group_by(platform) |>
  summarise(pct_hashtag = round(100 * mean(has_hashtag), 1),
            median_with = median(likeCount[has_hashtag]),
            median_without = median(likeCount[!has_hashtag]),
            .groups = "drop") |>
  kable(caption = "Hashtag use and median likes by platform")

Hashtag use and median likes by platform

platform	pct_hashtag	median_with	median_without
Facebook	36.2	76	54.0
Instagram	55.5	154	95.0
LinkedIn	30.4	121	149.0
TikTok	85.2	7	11.5
Twitter/X	26.4	2	1.0
YouTube	0.0	NA	0.0

# does hashtag use relate to engagement? two groups, so a t-test /
# Mann-Whitney U comparison
posts_log |> mutate(has_hashtag = factor(has_hashtag)) |> levene_test(log_likes ~ has_hashtag)

## # A tibble: 1 × 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1     1  8999      70.7 4.74e-17

t.test(log_likes ~ has_hashtag, data = posts_log)        # Welch's t-test (parametric)

## 
##  Welch Two Sample t-test
## 
## data:  log_likes by has_hashtag
## t = -20.169, df = 4566.4, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group FALSE and group TRUE is not equal to 0
## 95 percent confidence interval:
##  -0.4176917 -0.3436825
## sample estimates:
## mean in group FALSE  mean in group TRUE 
##           0.7013558           1.0820429

wilcox.test(log_likes ~ has_hashtag, data = posts_log)   # Mann-Whitney U (non-parametric)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  log_likes by has_hashtag
## W = 6043194, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

# median likes per topic, to show direction
posts_with_topic |>
  group_by(topic) |>
  summarise(posts = n(), median_likes = median(likeCount), .groups = "drop") |>
  kable(caption = "Median likes by topic")

Median likes by topic

topic	posts	median_likes
1	1087	1
2	1471	5
3	1484	8
4	1535	7
5	1387	4
6	1803	1

# likes across the six topics - Kruskal-Wallis with Dunn's test as post-hoc
topic_test <- kruskal.test(likeCount ~ topic, data = posts_with_topic)
topic_test

## 
##  Kruskal-Wallis rank sum test
## 
## data:  likeCount by topic
## Kruskal-Wallis chi-squared = 664.65, df = 5, p-value < 2.2e-16

dunn_test(posts_with_topic, likeCount ~ topic, p.adjust.method = "bonferroni") |>
  filter(p.adj < 0.05) |>
  select(group1, group2, statistic, p.adj, p.adj.signif)

## # A tibble: 13 × 5
##    group1 group2 statistic    p.adj p.adj.signif
##    <chr>  <chr>      <dbl>    <dbl> <chr>       
##  1 1      2          13.2  1.45e-38 ****        
##  2 1      3          19.0  3.96e-79 ****        
##  3 1      4          17.4  8.22e-67 ****        
##  4 1      5          12.2  8.41e-33 ****        
##  5 1      6           2.94 4.94e- 2 *           
##  6 2      3           6.25 6.17e- 9 ****        
##  7 2      4           4.47 1.19e- 4 ***         
##  8 2      6         -11.8  5.43e-31 ****        
##  9 3      5          -7.11 1.80e-11 ****        
## 10 3      6         -18.4  2.15e-74 ****        
## 11 4      5          -5.36 1.27e- 6 ****        
## 12 4      6         -16.6  5.54e-61 ****        
## 13 5      6         -10.6  3.51e-25 ****

# post length vs likes within each platform
posts |>
  group_by(platform) |>
  summarise(rho = round(cor(word_count, likeCount, method = "spearman"), 2),
            .groups = "drop") |>
  kable(caption = "Spearman rho: word count vs likes, by platform")

Spearman rho: word count vs likes, by platform

platform	rho
Facebook	0.10
Instagram	-0.18
LinkedIn	0.20
TikTok	0.10
Twitter/X	0.38
YouTube	0.15

# correlation matrix across the engagement metrics and post length,
# visualised as a correlogram
posts |>
  select(likeCount, commentCount, shareCount, videoplayCount, word_count) |>
  cor(method = "spearman", use = "pairwise.complete.obs") |>
  corrplot(type = "upper", tl.col ="black", addCoef.col ="black", number.cex = 0.9,
           col= colorRampPalette(c("#d73027", "white", "#1a9850")) (200),
           title = "Spearman correlations: engagement metrics and post length", mar = c (0,0,2,9))

# video plays vs. likes specifically
plays_likes_test <- cor.test(videos_with_plays$videoplayCount, videos_with_plays$likeCount,
                              method = "spearman", exact = FALSE)
plays_likes_test

## 
##  Spearman's rank correlation rho
## 
## data:  videos_with_plays$videoplayCount and videos_with_plays$likeCount
## S = 31148559, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.7312871

#bivariate_summary
# summary table of the headline result from each test
tibble(
  analysis = c("Likes by platform (Kruskal-Wallis)",
               "Likes by media type, high vs. low engagement (chi-squared)",
               "Likes by hashtag use (Mann-Whitney U)",
               "Likes by topic (Kruskal-Wallis)",
               "Post length vs. likes (Spearman)",
               "Video plays vs. likes (Spearman)"),
  statistic = c(paste0("H = ", round(unname(platform_test$statistic), 2)),
                paste0("X-squared = ", round(unname(chisq.test(chi_table1)$statistic), 2)),
                paste0("W = ", format(unname(wilcox.test(log_likes ~ has_hashtag, data = posts_log)$statistic), big.mark = ",")),
                paste0("H = ", round(unname(topic_test$statistic), 2)),
                paste0("rho = ", round(unname(cor.test(posts$word_count, posts$likeCount, method = "spearman", exact = FALSE)$estimate), 3)),
                paste0("rho = ", round(unname(plays_likes_test$estimate), 3))),
  p_value = c(platform_test$p.value,
              chisq.test(chi_table1)$p.value,
              wilcox.test(log_likes ~ has_hashtag, data = posts_log)$p.value,
              topic_test$p.value,
              cor.test(posts$word_count, posts$likeCount, method = "spearman", exact = FALSE)$p.value,
              plays_likes_test$p.value)
) |>
  mutate(p_value = format.pval(p_value, digits = 3, eps = 0.001)) |>
  kable(caption = "Bivariate test results - summary")

Bivariate test results - summary

analysis	statistic	p_value
Likes by platform (Kruskal-Wallis)	H = 4728.39	<0.001
Likes by media type, high vs. low engagement (chi-squared)	X-squared = 36.02	<0.001
Likes by hashtag use (Mann-Whitney U)	W = 6,043,194	<0.001
Likes by topic (Kruskal-Wallis)	H = 664.65	<0.001
Post length vs. likes (Spearman)	rho = 0.464	<0.001
Video plays vs. likes (Spearman)	rho = 0.731	<0.001

Commentary on engagement insights and recommendations

There are many more outliers than one would like - some posts have huge engagement compared to the rest, which means that medians will be used to describe typical engagement. Also, whenever the results from the Kruskal-Wallis test and Mann-Whitney test differ from the results of the parametric test, they are treated as primary.

Platform. The platforms with the highest median number of likes are LinkedIn (median = 141) and Instagram (median = 132), then Facebook (median = 61). According to both ANOVA (F = 3257, p < .001) and the Kruskal-Wallis test (H = 4728.4, p < .001), there are significant differences between the platforms. Post-hoc tests (Tukey and Dunn) show every pair differs significantly except Instagram vs LinkedIn and Twitter/X vs YouTube. However, it should be kept in mind that it is not an engagement-rate ranking - there are no measures of followers/following, reach, paid promotion and time lived in the data.

Media type. Image and video posts on Facebook, Instagram and LinkedIn were found to be more often than expected above the median number of likes for that platform compared to posts with not recorded media type (chi-square = 36.0, df = 3, p < .001; standardized residuals: Image +2.9, Video +3.2, Not recorded – 5.9). Sidecar posts had no significant difference. “Not recorded” refers to the media type being unknown rather than a text post, thus it gives rise to further experiments on image and video posts.

Hashtags. For all platforms combined, there are more likes for posts with hashtags compared to those without (Mann-Whitney U, W = 6,043,194, p < .001). However, hashtags vary across platforms, ranging from 0% on YouTube and 26–36% on Twitter/X, LinkedIn, and Facebook, to 56% on Instagram and 85% on TikTok, therefore some of the difference observed above might be attributable to the platform rather than hashtags. Looking at individual platforms, the trend is inconsistent: posts with hashtags get more likes on Facebook (median: 76 vs. 54), Instagram (154 vs. 95) and almost on Twitter/X (2 vs. 1), while less on LinkedIn (121 vs. 149) and TikTok (7 vs. 11.5). Medians above are descriptive only and have not been tested; therefore, hashtags cannot be considered an absolute engagement predictor.

Post length. When we consider all the platforms together, we observe a moderate positive correlation between caption length and likes (Spearman rho = 0.46, p < .001). In each platform separately, the correlation varies from -0.18 (Instagram) to 0.38 (Twitter/X): it is low but positive for Facebook (0.10), TikTok (0.10), YouTube (0.15), and LinkedIn (0.20), while it is negative for Instagram. There is variability in the strength and nature of the relationship, which indicates that platform context may be significant to understand the relationship between caption length and likes.

Video plays. Of 886 posts with plays information, there is a strong correlation between video plays and likes (Spearman rho = 0.73, p < .001). Again, correlation doesn’t imply causation here.

Topic. Likes vary greatly based on the topic (H for Kruskal-Wallis = 664.65, p < .001). The topic that has the greatest median likes is Topic 3 (community and schools) with a median value of 8, followed by Topic 4 (7, awards and recognition), Topic 2 (5) and Topic 5 (4). The topics having the fewest median likes are Topic 1 (business and pricing) and Topic 6 (team and events) at 1. In Dunn’s test (Bonferroni), most topic pairs showed significant differences; the non-significant pairwise differences were between Topics 2 and 5 and Topics 3 and 4. This analysis did not take into account the platforms and could also include specific platforms, so the next step would be comparing the topics in one platform over time.

Recommendations. In future data collection, Fonterra should include reach/impressions, follower count at the time of posting, paid/organic status and more complete share and view counts so that engagement rates can be calculated and platforms can be compared more fairly. Based on the current data, the results suggest several areas for further testing: (1) develop platform-specific strategies rather than relying on one approach, because median likes differ substantially across platforms; (2) test image and video content on Facebook, Instagram and LinkedIn, where these media types were more frequently associated with above-median likes; (3) test hashtag use separately by platform, particularly on Facebook and Instagram, while recognising that the observed differences do not establish that hashtags cause higher engagement; (4) avoid relying on caption length alone as an engagement strategy because the relationship varies substantially across platforms; (5) investigate community/school and awards/recognition content further, as these topics had higher median likes in this dataset; and (6) use rate-based engagement measures once reach and follower data are available.

Declaration for Generative AI use

• AI tool: (Copilot).
• Purpose of using the tool: I used AI to outline the analysis, to write and debug R code, and to ensure my commentary corresponded to the knitted results.
• Amount of AI usage: Copilot helped me to decide on which decision is important for the analysis and gave me the code as per my instruction and draft commentary. However, all the commentary and final decision were mine. I ran, reviewed, edited and made the final analytical decisions myself.
• Methods used to ensure accuracy:I knitted the document to make sure it works without any errors, compared all numbers in the commentary against the knitted results, read samples of posts on each topic to verify that the topic names were correct, and understood why each statistical test was used.