What happend in web tracking around 2022?

Introduction

Hello dear reader!

By this time, you probably know me reasonably well from all the CVs I have sent you. It seems like those weren’t that interesting to you, which I totally understand. But at the same time, CV is not really able to fully communicate how I actually work and whether the things that I’m able to create could add value to DuckDuckGo. That’s why I have decided to show you a piece of my work by collecting and analyzing your data from the tracker-radar project. I know that this might be a little bit weird and annoying. But, I see working at DuckDuckGo as a great opportunity - far too great to give up right away. With that being said, if you decide to reject my application again, after seeing a piece of my work, I will fully accept that and I won’t reach out again. Thank you!

– Stanislav “Stano” Smatana

What happend in 2022? Analysis of tracker-radar data in time.

tldr; It seems like something interesting have happend some time between 2022 and 2023 in the world of website trackers. The nature of trackers have suddenly shifted and this shift seems to be connected with fingerprinting based on viewport characteristics. Do you know what could have happend? If you want to get right into results, jump to the section Analyzing Suspected Tracker Patterns in Time. If you want to see how I went over analyzing this data step by step - continue reading.

I’ve decided to download and analyze website tracking data collected by DuckDuckGo’s tracker-radar with following goals:

• Characterize potential trackers based on combinations of API calls they are performing.
• Understand how character of trackers differs across big players and how it evolved in time.

All of the data I’ve used in the analysis was extracted from domain json files availiable across all releases of . You can see the shape of the raw input data below. Each records corresponds to one type of API call in one resource of one domain, meaning that the overall grain of the table is release_tag (date), region, domain, resource, api call. The data contains following columns:

• date - date of data collection, extracted from release_tag
• region - region in which this record was collected (corresponds to subfolders of domain folder)
• rule - regex identifying the resource
• api - api function
  
• calls - number of calls to the api function
• fingerprinting - fingerprinting score derived in tracker-radar for the resource

library(dplyr)
library(topicmodels)

api_calls_in_time <- read.csv("data/api_calls_by_release.csv") %>%
  filter(type == "Script") 

#converting release tags to dates
release_tag_to_date <- data.frame(release_tag = unique(api_calls_in_time$release_tag)) %>%
  separate(release_tag, into = c("year", "month", "day"), sep = "\\.", fill = "right",remove = F) %>%
  mutate(
    month = sprintf("%02d", as.integer(month)),
    day   = ifelse(is.na(day), "01", sprintf("%02d", as.integer(day))),
    date  = as.Date(paste(year, month, day, sep = "-"))
  ) %>%
  select(release_tag,date)

api_calls_in_time <- api_calls_in_time %>%
  inner_join(release_tag_to_date) %>%
  select(date,region,domain,rule,api,calls,fingerprinting)
  
head(api_calls_in_time)

Python scripts I’ve used to scrape tracker-radar together with the source code of this notebook and other artifacts are availiable here TODO.

NOTE: Python is my primary language, but I love to use R for exploratory data analysis, so that’s what I’m using here :)

AI Use During this Analysis

During this process I’ve used AI (GitHub Copilot) at times to complement my primary skills. I have used AI exactly in this way:

• I have used AI to generate python scripts for scraping of , because I didn’t want to waste time implementing scraping and to focus on analysis.
• I did not use AI to conduct the analysis of this data and create this notebook, because that’s the part that I enjoy, where I excel, and where I do a better job.
• I have used AI to help interpret LDA topics (details in the next section), because my knowledge of web trackers is very limited :)

Characterizing API Calls - Latent Dirichlet Allocation

Fingerprinting is connected to certain Javascript API calls that can help create precise identification of a device or a user. My assumption is that there are certain types of trackers and that each tracker type is associated with certain API calls. In order to see what these suspected tracker types might be, I need first to identify API calls that tend to occur together - API call groups.

Nearly any clustering approach could be used to create some groups of API calls, but I have decided to use Latent Dirichlet Allocation (LDA, not to be confused with Linear Discriminant Analysis). While this model traditionally comes from the field of Natural Language Processing, its assumptions exactly match this setting - sparse count data (counts of api calls with only few appearing in each resource). Don’t be scared by the complicated name, it’s actually a fairly simple unsupervised generative model. Let me explain how it specifically works for this analysis.

Every statistical model assumes certain view of the world and then helps us derive new information based on this view. In this setting, LDA assumes that certain API calls tend to occur together. These co-ocuring calls can be grouped into larger call groups. Each call can be a member of multiple groups and this membership is soft. This means that each call is associated with each group (topic in LDA language) with certain probability between 0 and 1.

With groups established, script resources (documents in LDA language) are viewed as different mixes of groups. Each mix is characterized by resource to group associations between 0-1 that add up to 1. Groups, associations between resources and groups and associations between groups and calls are derived by LDA from data without supervision. The only required parameter is the number of groups (topics). I will describe selection of the number of groups later.

The diagram below illustrates how LDA models an example of 2 script resources. Each script is associated with call groups 1-4 and its group mix dictates the actual calls that appear in each script resource. For instance, example1.js is 70% group1 and 70% group2. Therefore, its API calls come mostly from group2 and then also from group1. This is a toy example, in real LDA, each document is usually associated with each topic to an extent.

What does LDA enable me to do?

• Reduce dimensionality in an interpretable way. Instead of characterizing resources by all api call types, I can characterize them by whole call groups.
• Understand which API calls tend to co-occur and understand what what these groups of co-occuring calls might represent.

In order to derive call groups I will use only data from the most recent release. Afterwards, I’ll apply the model to the full data spanning all tracker-radar releases. I do this because I assume that today we have the broadest spectrum of tracking techniques available and the past was equally or less diverse.

api_calls_recent <- 
  api_calls_in_time %>%
  filter(date == max(api_calls_in_time$date))

head(api_calls_recent)

One last thing I need to do in order to apply LDA is to determine the number of call groups (topics) - k. There are many methods of doing that. I have decided to use one based on coisne similarity (Cao et. al. 2009), which aims to create groups that are as distinct from one another as possible. This is helpful because LDA sometimes tends to create groups (topics) that are heavily overlapping. Presence of multiple overlapping groups that correspond roughly to the same concept creates noise that makes interpretation harder.

More formally, k, for which the model has groups with smallest average pairwise cosine similarity is selected as the winning k.

cosine_sim <- function(M) {
  # row norms
  norms <- sqrt(rowSums(M * M))
  
  # denominator via outer
  den <- outer(norms, norms, "*")
  
  # dot products
  dot <- M %*% t(M)
  
  # cosine similarity
  cos_sim <- dot / den
  cos_sim
}

avg_cos_sim <- function(M){
  csim <- cosine_sim(M)
  mean(csim[upper.tri(csim)])  
}

call_counts <- api_calls_recent %>%
  pivot_wider(names_from=api,values_from = calls,values_fill = 0) 

call_counts_mat <- call_counts %>% 
  select(-c(region,domain,rule,fingerprinting,date)) %>% 
  as.matrix()

#This takes a lot of time, just load the rds if winner was already calculated
if(file.exists("data/winning_lda.Rds")){
  winning_lda = readRDS("data/winning_lda.Rds")
}else
{
  csims <- rep(0,times=19)
  for(k in 2:20){
    cat(paste("Processing LDA ",k,"\n"))
    topic_model <- LDA(call_counts_mat,k)    
    csims[k-1] <- avg_cos_sim(topic_model@beta)
  }

  
  selected_k <- (2:20)[which.min(csims)]
  winning_lda = LDA(call_counts_mat,selected_k)
  
  
}

selected_k = ncol(winning_lda@gamma)

LDA with 4 topics was selected as the model with the lowest average cosine similarity between API call groups (topics). Next, I want to characterize each group (topic) and for that I’ll identify calls that are most characteristic for each group. This can be done by looking into the beta matrix that has 4 rows, one for each group, and as many columns as there are distinct call types. Each cell corresponds to probability (in log scale) of corresponding call being associated with a corresponding call group. Small excerpt of the matrix below:

assoc_mat <- winning_lda@beta[,1:3]
colnames(assoc_mat) <- winning_lda@terms[1:3]
rownames(assoc_mat) <- paste("Group",1:4)
assoc_mat

        HTMLCanvasElement.prototype.toDataURL Navigator.prototype.hardwareConcurrency Navigator.prototype.plugins
Group 1                            -46.927331                              -39.063618                  -27.901858
Group 2                             -9.249582                              -25.311132                  -20.326818
Group 3                            -31.002229                              -12.835772                  -31.452065
Group 4                             -5.318555                               -3.823308                   -3.639734

However, it is important to note that if a call has high probability of appearing in a given group, it doesn’t automatically mean that it is very characteristic of this group. It might be the case that it is a very common call that has equally high probability of appearing in each group.

Luckily, it is possible to use Bayes rule to calculate the probability of call being characteristic of a given group. This is done simply by normalizing the columns of the beta matrix into 0-1 range. I have used logsumexp and subtraction before applying exp. The same result could be achieved by exponentiating first and then normalizing using division by sum of probabilities. However, calculations in log scale are numerically more stable, so it makes sense to stay in the log scale as long as possible.

Next, I identify the most characteristic API calls for each group - API calls that have probability of being characteristic of a given group higher than 80%.

characteristic_terms <- function(lda_model){
  
  sums <- apply(lda_model@beta,2,logSumExp)
  p_distinctive <- exp(sweep(lda_model@beta,2,sums))  
  
  n_topics <- nrow(lda_model@beta)
  results <- list()
  
  for(i in 1:n_topics){
    ordering <- order(lda_model@beta[i,],decreasing = T)
    ordered_terms <- lda_model@terms[ordering]
    ordered_selection <- (distinctive[i,] > 0.8)[ordering]
    selected_ordered_terms <- ordered_terms[ordered_selection]
    ordered_ps <- (exp(lda_model@beta[i,])[ordering])[ordered_selection]      
    names(ordered_ps) <- selected_ordered_terms
    
    results[[i]] <- ordered_ps
  }
  
  results

}

characteristic_calls <- characteristic_terms(winning_lda)

Because my knowledge of web tracking methods is fairly limited, I have consulted AI at this stage. For each group, I have taken the characteristic API calls and asked AI (Copilot) these questions:

• “What these javascript function calls have in common?”
• “How would you characterize this group?”
• “What name would you give to this group of functions?”

I have made sure to ask these questions without mentioning web tracking. Interestingly, for each group, the AI have immediately identified that tracking/fingerprinting might be at play. Answers helped me establish names and characterization of these API call groups.

The series of barplots below shows characteristic calls for each call group. The length of the bar corresponds to how frequently is given call associated with a given group.

Identity & Persistance API Call Group

Functions focused on persistante storage and high-entropy identifiers.

Browser Personality API Call Group

Calls focused on characterizing user’s browser. This is a broad group - limiting to top 15 most common calls.

Rendering & Viewing Surface API Call Group

Mostly properties of the viewport, rendering quirks and motion sensor.

Capability & Hardware API Call Group

Mostly functions that can be used to get broad information on hardware and device capabilities. Big group, showing only top 15.

API Call Groups and Suspected Tracking

Application of LDA above have helped to identify certain groups of calls that tend to co-occur. This does not mean, however, that all of them are indication of fingerprinting. In order to see which groups or group combinations might be associated with fingerprinting, I’m using logistic regression on resource-group associations, with outcome (predicted) variable being the fingerprinting score calculated by tracker-radar. I want to see which call groups, or their combinations are associated with high fingerprinting score.

Of course, the fingerprinting score itself is just a heuristic and might not be always correct. But this is ok for initial exploratory analysis to discover interesting hypotheses. To make things easier for interpretation, I have made two simplifications:

• I have binarized associations between call groups and script resources. Call group is associated with a given script resource, if its association probability with given group is greater than what would be expected in 70% of cases.
• Instead of predicting fingerprinting score of 1 - 3, I’m predicting binary outcome of fingerprinting score being 1 vs 2 or 3.

Remember that I’m not using classifier here to accurately predict the fingerprinting score. The main purpose is to see which call groups or call group combinations might be associated with high fingerprinting scores and therefore likely associated with trackers/fingerprinting. I’m also not doing the more common train-test-validation split, but using information criterion (AIC) for model selection. I’m not interested in how exactly this model would work in production - I’m just interested in choosing better model.

You can see the excerpt of data used for training below:

assoc_thresholds <- colQuantiles(lda_winner@gamma,probs = 0.7)
binary_activations <- lda_winner@gamma > assoc_thresholds
binary_activations_df <- as.data.frame(activations)
binary_activations_df$outcome <- call_counts$fingerprinting > 1

head(binary_activations_df)

I have tried 2 versions of logistic regression:

• Regressing on each variable, but not on their combinations
• Regression on each variable and each pairwise variable combination (pair interactions)

The code below fits these models and outputs their AIC score (lower is better).

pred_model_base <- glm(outcome ~ V1 + V2 + V3 + V4,binary_activations_df,family = "binomial")
pred_model_pairs <- glm(outcome ~ V1 + V2 + V3 + V4 + V1*V2 + V1*V3 + V1*V4 + V2*V3 + V2*V4 + V3*V4,binary_activations_df,family = "binomial")
data.frame(model=c("Baseline Model","Interaction Model"),AIC = c(AIC(pred_model_base),
AIC(pred_model_pairs)))

Interaction model is a clear winner. The difference in AIC value might seem small at the first sight. However, AIC is defined on a logarithmic scale, so this difference is actually very significant. Next, I have generated all possible combinations of activations to see which will be predicted by the model to have more than 50% probability of tracking behavior.

The table below shows the call group combinations associated with tracking/fingerprinting behavior (p>50%) ordered from most probable to the least probable. Acronym for each combination was created by taking first letters of each call group taht it consists of.

all_combs <- expand.grid(V1=c(T,F),V2=c(T,F),V3=c(T,F), V4 = c(T,F))
all_combs$Probability <- predict(pred_model_pairs,all_combs,type="response")


all_combs %>%
  filter(Probability > 0.5) %>%
  arrange(desc(Probability)) %>%
  mutate(
    Acronym = paste0(ifelse(V1,"I",""),ifelse(V2,"BP",""),ifelse(V3,"R",""),ifelse(V4,"C","")),
    V1 = ifelse(V1,"Identity&Persistance",""),
    V2 = ifelse(V2,"Browser Personality",""),
    V3 = ifelse(V3,"Rendering&Viewport",""),
    V4 = ifelse(V4,"Capability","")
  ) %>%
  select(Acronym,V1,V2,V3,V4,Probability)

Analyzing Suspected Tracker Patterns in Time

Having established the API call groups and their combinations associated with fingerprinting, I can now use these combinations to see how the character of fingerprinting evolved in time. First, I’ll investigate overall development and then I’ll investigate the biggest players in tracking.

In order to do that, I had to do some data processing focused on annotation of each script resource by the call group combination it belongs to.

call_counts_in_time <-api_calls_in_time %>%
  pivot_wider(names_from=api,values_from = calls,values_fill = 0) 

call_counts_mat_in_time <- call_counts_in_time %>% 
  select(-c(date,region,domain,rule,fingerprinting)) %>% 
  as.matrix()

#Order columns in the same way as in LDA training
call_counts_mat_in_time <- call_counts_mat_in_time[,colnames(call_counts_mat)]

non_zero_row <- rowSums(call_counts_mat_in_time) > 0
call_counts_mat_in_time <- call_counts_mat_in_time[non_zero_row,]
call_counts_in_time_nonzero <- call_counts_in_time[non_zero_row,]

topic_posteriors <- posterior(lda_winner,call_counts_mat_in_time)
activations_in_time <- topic_posteriors$topics > assoc_thresholds
colnames(activations_in_time) <- c("V1","V2","V3","V4")

tracker_types_in_time <- activations_in_time %>% as.data.frame %>%
  mutate(
    tracker_type = case_when(
      V1 & V2 & V3 & V4 ~ "IBPRC",   # Identity + Browser Personality + Rendering + Capability
      V2 & V3 & V4      ~ "BPRC",    # Browser Personality + Rendering + Capability
      V2 & V4           ~ "BPC",     # Browser Personality + Capability
      V3 & V4           ~ "RC",      # Rendering + Capability
      V4                ~ "C",       # Capability only
      V1 & V2 & V3      ~ "IBPR",    # Identity + Browser Personality + Rendering
      V1 & V2 & V4      ~ "IBPC",    # Identity + Browser Personality + Capability
      V1 & V3 & V4      ~ "IRC",     # Identity + Rendering + Capability
      TRUE              ~ NA_character_
    )
  ) 

tracker_types_in_time$domain <- call_counts_in_time_nonzero$domain
tracker_types_in_time$date <- call_counts_in_time_nonzero$date

The plot below shows how the number of scripts associated with each pattern evolved over time. There are a couple of things to note here:

• Out of the 10 patterns suggesting tracking according to the logistic regression model, only 5 are actually really occurring in the wild. This is expected behavior - statistical models can be used for speculation over fictional scenarios. I have limited further analysis to these real patterns:
  • BPC - Browser Personality and Capability
  • BPRC - Browser Personality, Rendering and Capability
  • C - Capability
  • IBPR - Identity, Browser Personality and Rendering&Viewport
• The plot clearly shows that the predominant pattern is Capability-only followed by the combination of Browser Personaity and Capability.
• However, at some time between 2022 and 2023 there’s a sudden uptick in two patterns with functions known to be involved in viewport fingerprinting - IBPR and especially RC.

The plot below shows the same evolution of suspected tracker patterns in both absolute and relative terms (percentages) broken down for each major player. It is apparent that the same pattern seems to hold with nearly all major players. Especially Google, Akamai and Adobe seem to have sudden uptick in IBPR and RC, while Microsoft seems to always have had a significant part of its resources aligned with the RC pattern.

Finally, I’ve decided to zoom in into the grain of individual calls. I wanted to answer the question - how do the proportions of different API calls differ pre and post 2022? In order to identify changes that are significant, I’m using chi.square test over 2x2 contingency matrices. This means that for every call I create a matrix of the form

X	Calls_of_function_f	Other_calls
Before 22	n_f_pre22	n_other_pre22
After 22	n_f_post22	n_other_post22

And then apply chis.square test to derive significance p-values. I’ve chosen the significance threshold 0.01. However, because I’m performing a great number of these tests, the significance threshold needs to be divided by the number of tests so that results have expected error rate (this is called Bonferroni correction). Commentary of results is part of the next section - Conclusion.

period_api_calls <- api_calls_in_time %>% 
  select(-c(region,domain,rule)) %>%
  mutate(pre_2022 = date < as.Date("2022-1-1")) %>%
  group_by(pre_2022,api) %>%
  summarize(calls = sum(calls))

period_api_calls_wide <- period_api_calls %>%
  group_by(pre_2022) %>%
  mutate(period_total_calls = sum(calls)) %>%
  ungroup() %>%
  pivot_wider(names_from = pre_2022,values_from = c(calls,period_total_calls),values_fill = 0) %>%
  mutate(period_total_calls_TRUE = ifelse(period_total_calls_TRUE == 0,period_total_calls_TRUE[period_total_calls_TRUE > 0][1],period_total_calls_TRUE))


chi_square_calls <- function(df){
  calls_pre22 <- df$calls_TRUE
  non_calls_pre22 <- df$period_total_calls_TRUE - df$calls_TRUE
  calls_post22 <- df$calls_FALSE
  non_calls_post22 <- df$period_total_calls_FALSE - df$calls_FALSE
  
  all_data_mat <- cbind(calls_pre22,non_calls_pre22,calls_post22,non_calls_post22)
  apply(all_data_mat,1,function(row){
    chisq.test(rbind(
      c(row["calls_pre22"],row["non_calls_pre22"]),
      c(row["calls_post22"],row["non_calls_post22"])
    ))$p.value
  })
}

period_api_calls_wide$significant <- chi_square_calls(period_api_calls_wide) < 7.874016e-05 #bonferoni correction
period_api_calls_wide %>% 
  filter(significant) %>%
  select(api,calls_TRUE,calls_FALSE) %>%
  arrange(calls_TRUE) %>%
  rename(calls_before_2022 = calls_TRUE,calls_after_2022 = calls_FALSE)

Conclusion

It does seem like on top of the shift in call groups, there is a significant difference in proportions before 2022 and after 2022 for pretty much every individual API call. What’s more important some heavily used API calls are completely absent before 2022. This leads me to a couple of hypothesis explaining the sudden shift of tracking patterns in 2022:

• There was a true change of how companies are doing fingerprinting motivated by new approaches or change in (legislative?) environment.
• Some API calls useful for fingerprinting were not available in JS before 2022. Their addition motivated companies to include them in their tracking approaches.
• The way tracker-radar collects data have changed and there’s no real change in behavior of organizations that do fingerprinting/tracking.
• Or any combination of above.

I have a question for you, dear reader, what do you think is the most likely explanation? Feel free to share any ideas you have on my linkedin or via my e-mail stanislavsmatana@gmail.com.