This does indeed appear to be a bug in Plotly - this can be submitted as a bug report to the Plotly team.
It is worth noting that modifying boxpoints = "outliers" to boxpoints = "suspectedoutliers" produces markers with a different color so suspectedoutliers behaves as expected. However, you can't use suspectedoutliers in place of outliers as suspected outliers are only a subset of all outliers.
You can achieve the desired behavior by plotting the outliers manually. To do this, you would still set boxpoints=outliers, but then plot the outliers as individual scatter points with the desired color over the outliers generated by Plotly.
This is a bit intensive because this requires a rewrite of the algorithm to determine outliers exactly as the Plotly library performs this calculation. And unfortunately, you cannot extract Q1, Q3 or other statistics from go.Box or from Plotly in any way as these computations are performed by the Javascript under the hood when the figure renders.
The first thing to note is that calculating Q1 and Q3 differs between different Python libraries: Plotly outlines their methods in the documentation, explaining that they use Method #10 in this short paper to calculate percentiles.
In Python, the function to calculate percentiles using Method #10 (linear interpolation) looks like this:
## calculate quartiles as outlined in the plotly documentation
def get_percentile(data, p):
data.sort()
n = len(data)
x = n*p + 0.5
x1, x2 = floor(x), ceil(x)
y1, y2 = data[x1-1], data[x2-1] # account for zero-indexing
return y1 + ((x - x1) / (x2 - x1))*(y2 - y1)
Now to extract outliers from a data set, you subset the data: anything below (Q1 - 1.5 * IQR) or above (Q3 + 1.5 * IQR) where IQR = Q3 - Q1 is considered an outlier.
Putting this all together:
from math import floor, ceil
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from matplotlib.colors import LinearSegmentedColormap, to_hex
df_plot = pd.read_csv('https://raw.githubusercontent.com/mwaskom/seaborn-data/master/iris.csv')
cat_var = "species"
num_var = "petal_length"
lvls = df_plot[cat_var].unique()
n_levels = len(lvls)
cmap = LinearSegmentedColormap.from_list("my_palette", ["#111539", "#97A1D9"])
my_palette = [to_hex(j) for j in [cmap(i/n_levels) for i in np.array(range(n_levels))]]
## calculate quartiles as outlined in the plotly documentation
def get_percentile(data, p):
data.sort()
n = len(data)
x = n*p + 0.5
x1, x2 = floor(n*p), ceil(n*p)
y1, y2 = data[x1-1], data[x2-1] # account for zero-indexing
return y1 + ((x - x1) / (x2 - x1))*(y2 - y1)
def get_fences(data):
q1, q3 = get_percentile(data, 0.25), get_percentile(data, 0.75)
iqr = q3-q1
return (q1 - (1.5*iqr), q3 + (1.5*iqr))
boxes = []
for l in range(n_levels):
data = df_plot.loc[df_plot.loc[:, cat_var] == lvls[l], num_var].values
outliers = data[(data < get_fences(data)[0]) | (data > get_fences(data)[1])]
print(outliers)
boxes += [
go.Box(
name = lvls[l],
y = data,
width = 0.4,
boxpoints = "outliers",
marker = {
"outliercolor": "red", ### there may be a plotly.go bug here
"color": my_palette[l],
"size": 30,
"opacity": 0.5
}
),
go.Scatter(
x = [lvls[l]]*len("outliers"),
y = outliers,
mode = 'markers',
marker=dict(color="red", size=28, opacity=0.5)
)
]
fig = go.Figure(data = boxes)
fig.update_layout(
font = dict(
size = 18
),
showlegend = False,
plot_bgcolor = "white",
hoverlabel = dict(
font_size = 18,
font_family = "Rockwell"
)
)
fig.show()

As a way to check our work, you will notice that the slightly smaller manually added outliers match the outliers determined by Plotly. (You can make the manually added outliers larger to obscure the Plotly generated outliers that aren't the desired color)