32:T41c,grid = Plot.plot({ width: 600, padding: 0, legend: true, grid: true, x: { label: "Episode", axis: "top" }, y: { label: "Season" }, color: { label: "Rating", type: "linear", scheme: "RdYlGn", legend: true, domain: [0, 10], unknown: "#ddd" }, marks: [ Plot.cell(cohortAllCohortAverages, { x: "episode", y: "season", fill: "cohortAverageRating", inset: 0.5, tip: true, channels: { Title: (d) => (d.title ? d.title : "???") } }), Plot.text( cohortAllCohortAverages, colorContrast({ x: "episode", y: "season", text: (d) => d.cohortAverageRating?.toFixed(1), fill: "cohortAverageRating" }) ), Plot.ruleY([cohortAllCohortAverages[goldenAgeFinalEpisodeIndex]], { x: "episode", y: "season", dy: 10, stroke: "red" }), Plot.dot([cohortAllCohortAverages[goldenAgeFinalEpisodeIndex]], { x: "episode", y: "season", dy: 10, fill: "red" }) ] })33:T4b3,mainPlot = Plot.plot({ width: Math.min(width, 640), color: { label: "Rating", type: "linear", scheme: "RdYlGn", legend: true, domain: [0, 10], unknown: "#ddd" }, r: { range: [0, 10] }, x: { domain: [0, 0.3], label: "Proportion of words spoken by Lisa (%) →" }, marginRight: 83, fx: { domain: ["males", "all", "females"] }, fy: { domain: ["all", "aged_under_18", "aged_18_29", "aged_30_44", "aged_45_plus"] }, marks: [ Plot.frame(), // Plot.linearRegressionY(cohortAverages, { // x: "lisaProportion", // y: "cohortAverageRating", // fx: "gender", // fy: "age" // }), Plot.linearRegressionY(exploded, { // ci: 1 - 0.05 ** 10, (TODO: ASK HUW) x: "lisaProportion", y: "rating", fx: "gender", fy: "age" }), Plot.dot(cohortAverages, { x: "lisaProportion", y: "cohortAverageRating", fx: "gender", fy: "age", opacity: 0.2, r: "cohortTotalVotes", fill: "cohortAverageRating", stroke: "black", tip: true, channels: { S_E_: "seNum", episode: "ep", cohort: "cohort" } }) ] })34:T533,```r # Install rjson if not already installed if (!require("rjson", quietly = TRUE)) { install.packages("rjson", repos = "https://cloud.r-project.org") library(rjson) } # Read the exploded data from JSON using readLines for better performance cat("Reading JSON file...\n") json_lines <- readLines("exploded.json") json_str <- paste(json_lines, collapse = "") cat("Parsing JSON...\n") json_data <- fromJSON(json_str) cat("JSON loaded:", length(json_data), "entries\n\n") # Convert list to data frame, handling NULLs lisaProp <- unlist(lapply(json_data, function(x) if(is.null(x$lisaProportion)) NA else x$lisaProportion)) ratings <- unlist(lapply(json_data, function(x) if(is.null(x$rating)) NA else x$rating)) ages <- unlist(lapply(json_data, function(x) if(is.null(x$age)) NA else x$age)) genders <- unlist(lapply(json_data, function(x) if(is.null(x$gender)) NA else x$gender)) df <- data.frame( lisaProportion = lisaProp, rating = ratings, age = ages, gender = genders, stringsAsFactors = FALSE ) # Remove any NA rows df <- df[complete.cases(df), ] # Convert to factors for modeling df$rating <- as.numeric(df$rating) df$lisaProportion <- as.numeric(df$lisaProportion) df$age <- as.factor(df$age) df$gender <- as.factor(df$gender) model <- lm(rating ~ lisaProportion * age * gender, data = df) anova(model) ```35:T3f82,function anova(exploded) { // Nested utility functions function mean(arr) { return arr.reduce((a, b) => a + b, 0) / arr.length; } function fCDF(x, df1, df2) { const z = (df1 * x) / (df1 * x + df2); return incompleteBeta(z, df1 / 2, df2 / 2); } function incompleteBeta(x, a, b) { if (x <= 0) return 0; if (x >= 1) return 1; const lbeta = logGamma(a) + logGamma(b) - logGamma(a + b); const front = Math.exp(Math.log(x) * a + Math.log(1 - x) * b - lbeta) / a; const f = betaCF(x, a, b); return front * f; } function betaCF(x, a, b, maxIter = 100) { const qab = a + b; const qap = a + 1; const qam = a - 1; let c = 1; let d = 1 - (qab * x) / qap; if (Math.abs(d) < 1e-30) d = 1e-30; d = 1 / d; let h = d; for (let m = 1; m <= maxIter; m++) { const m2 = 2 * m; let aa = (m * (b - m) * x) / ((qam + m2) * (a + m2)); d = 1 + aa * d; if (Math.abs(d) < 1e-30) d = 1e-30; c = 1 + aa / c; if (Math.abs(c) < 1e-30) c = 1e-30; d = 1 / d; h *= d * c; aa = (-(a + m) * (qab + m) * x) / ((a + m2) * (qap + m2)); d = 1 + aa * d; if (Math.abs(d) < 1e-30) d = 1e-30; c = 1 + aa / c; if (Math.abs(c) < 1e-30) c = 1e-30; d = 1 / d; const del = d * c; h *= del; if (Math.abs(del - 1) < 1e-10) break; } return h; } function logGamma(z) { const g = 7; const coef = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ]; if (z < 0.5) { return ( Math.log(Math.PI) - Math.log(Math.sin(Math.PI * z)) - logGamma(1 - z) ); } z -= 1; let x = coef[0]; for (let i = 1; i < g + 2; i++) { x += coef[i] / (z + i); } const t = z + g + 0.5; return ( 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x) ); } function formatPValue(p) { if (p < 2.2e-16) return "< 2.2e-16"; if (p < 0.001) { let exp = p.toExponential(3); // Convert e-7 to e-07 exp = exp.replace(/e-(\d)$/, "e-0$1"); return exp; } return p.toFixed(4); } function formatSigCode(p) { if (p < 0.001) return " ***"; if (p < 0.01) return " ** "; if (p < 0.05) return " * "; if (p < 0.1) return " . "; return " "; } function createDummyMatrix(variable, levels) { const dummies = []; for (let i = 1; i < levels.length; i++) { dummies.push(variable.map((v) => (v === levels[i] ? 1 : 0))); } return dummies; } function createInteractions(matrix1, matrix2) { const interactions = []; for (const col1 of matrix1) { for (const col2 of matrix2) { interactions.push(col1.map((v, i) => v * col2[i])); } } return interactions; } function gaussianElimination(A, b) { const n = A.length; const aug = A.map((row, i) => [...row, b[i]]); const eps = 1e-10; for (let i = 0; i < n; i++) { let maxRow = i; for (let k = i + 1; k < n; k++) { if (Math.abs(aug[k][i]) > Math.abs(aug[maxRow][i])) { maxRow = k; } } [aug[i], aug[maxRow]] = [aug[maxRow], aug[i]]; if (Math.abs(aug[i][i]) < eps) { continue; } for (let k = i + 1; k < n; k++) { const c = aug[k][i] / aug[i][i]; for (let j = i; j <= n; j++) { aug[k][j] -= c * aug[i][j]; } } } const x = Array(n).fill(0); for (let i = n - 1; i >= 0; i--) { if (Math.abs(aug[i][i]) < eps) { x[i] = 0; continue; } x[i] = aug[i][n]; for (let j = i + 1; j < n; j++) { x[i] -= aug[i][j] * x[j]; } x[i] /= aug[i][i]; } return x; } function fitModel(y, predictors) { const n = y.length; const p = predictors.length + 1; const X = Array(n) .fill(null) .map(() => Array(p).fill(0)); for (let i = 0; i < n; i++) { X[i][0] = 1; for (let j = 0; j < predictors.length; j++) { X[i][j + 1] = predictors[j][i]; } } const XtX = Array(p) .fill(null) .map(() => Array(p).fill(0)); for (let i = 0; i < p; i++) { for (let j = 0; j < p; j++) { for (let k = 0; k < n; k++) { XtX[i][j] += X[k][i] * X[k][j]; } } } const Xty = Array(p).fill(0); for (let i = 0; i < p; i++) { for (let k = 0; k < n; k++) { Xty[i] += X[k][i] * y[k]; } } const beta = gaussianElimination(XtX, Xty); const fitted = X.map((row) => row.reduce((sum, x, i) => sum + x * beta[i], 0) ); const rss = y.reduce((sum, yi, i) => sum + (yi - fitted[i]) ** 2, 0); return { beta, fitted, rss }; } // Main computation const rawData = exploded; const data = rawData.filter( (d) => d.rating != null && d.lisaProportion != null && d.age != null && d.gender != null && !isNaN(d.rating) && !isNaN(d.lisaProportion) ); const y = data.map((d) => d.rating); const lisaProp = data.map((d) => d.lisaProportion); const age = data.map((d) => d.age); const gender = data.map((d) => d.gender); const n = y.length; const yMean = mean(y); const TSS = y.reduce((sum, val) => sum + (val - yMean) ** 2, 0); const ageUnique = [...new Set(age)].sort(); const genderUnique = [...new Set(gender)].sort(); const ageMatrix = createDummyMatrix(age, ageUnique); const genderMatrix = createDummyMatrix(gender, genderUnique); const lisaAgeInteraction = createInteractions([lisaProp], ageMatrix); const lisaGenderInteraction = createInteractions([lisaProp], genderMatrix); const ageGenderInteraction = createInteractions(ageMatrix, genderMatrix); const lisaAgeGenderInteraction = createInteractions( [lisaProp], ageGenderInteraction ); let RSS0 = TSS; const model1 = fitModel(y, [lisaProp]); const RSS1 = model1.rss; const SS_lisaProp = RSS0 - RSS1; const model2 = fitModel(y, [lisaProp, ...ageMatrix]); const RSS2 = model2.rss; const SS_age = RSS1 - RSS2; const model3 = fitModel(y, [lisaProp, ...ageMatrix, ...genderMatrix]); const RSS3 = model3.rss; const SS_gender = RSS2 - RSS3; const model4 = fitModel(y, [ lisaProp, ...ageMatrix, ...genderMatrix, ...lisaAgeInteraction ]); const RSS4 = model4.rss; const SS_lisaAge = RSS3 - RSS4; const model5 = fitModel(y, [ lisaProp, ...ageMatrix, ...genderMatrix, ...lisaAgeInteraction, ...lisaGenderInteraction ]); const RSS5 = model5.rss; const SS_lisaGender = RSS4 - RSS5; const model6 = fitModel(y, [ lisaProp, ...ageMatrix, ...genderMatrix, ...lisaAgeInteraction, ...lisaGenderInteraction, ...ageGenderInteraction ]); const RSS6 = model6.rss; const SS_ageGender = RSS5 - RSS6; const model7 = fitModel(y, [ lisaProp, ...ageMatrix, ...genderMatrix, ...lisaAgeInteraction, ...lisaGenderInteraction, ...ageGenderInteraction, ...lisaAgeGenderInteraction ]); const RSS7 = model7.rss; const SS_lisaAgeGender = RSS6 - RSS7; const RSSResidual = RSS7; // Compute summary statistics for the full model const residuals = y.map((yi, i) => yi - model7.fitted[i]); const sortedResiduals = [...residuals].sort((a, b) => a - b); const resMin = sortedResiduals[0]; const res1Q = sortedResiduals[Math.floor(n * 0.25)]; const resMedian = sortedResiduals[Math.floor(n * 0.5)]; const res3Q = sortedResiduals[Math.floor(n * 0.75)]; const resMax = sortedResiduals[n - 1]; // Build coefficient names const coefNames = [ "(Intercept)", "lisaProportion", ...ageUnique.slice(1).map((a) => `age${a}`), ...genderUnique.slice(1).map((g) => `gender${g}`), ...ageUnique.slice(1).map((a) => `lisaProportion:age${a}`), ...genderUnique.slice(1).map((g) => `lisaProportion:gender${g}`), ...ageUnique .slice(1) .flatMap((a) => genderUnique.slice(1).map((g) => `age${a}:gender${g}`)), ...ageUnique .slice(1) .flatMap((a) => genderUnique.slice(1).map((g) => `lisaProportion:age${a}:gender${g}`) ) ]; // Compute standard errors for coefficients const nCoef = model7.beta.length; const dfModel = nCoef - 1; const df_residual_summary = n - nCoef; const sigma2 = RSSResidual / df_residual_summary; // Build X matrix for full model to compute (X'X)^-1 const XFull = Array(n) .fill(null) .map(() => Array(nCoef).fill(0)); const fullPredictors = [ lisaProp, ...ageMatrix, ...genderMatrix, ...lisaAgeInteraction, ...lisaGenderInteraction, ...ageGenderInteraction, ...lisaAgeGenderInteraction ]; for (let i = 0; i < n; i++) { XFull[i][0] = 1; for (let j = 0; j < fullPredictors.length; j++) { XFull[i][j + 1] = fullPredictors[j][i]; } } // Compute X'X const XtX_full = Array(nCoef) .fill(null) .map(() => Array(nCoef).fill(0)); for (let i = 0; i < nCoef; i++) { for (let j = 0; j < nCoef; j++) { for (let k = 0; k < n; k++) { XtX_full[i][j] += XFull[k][i] * XFull[k][j]; } } } // Compute (X'X)^-1 using Gaussian elimination on identity const XtX_inv = Array(nCoef) .fill(null) .map((_, i) => { const identity = Array(nCoef).fill(0); identity[i] = 1; return gaussianElimination(XtX_full, identity); }); // Standard errors are sqrt(diag(sigma2 * (X'X)^-1)) const stdErrors = XtX_inv.map((row, i) => Math.sqrt(sigma2 * row[i])); // t-values and p-values function tCDF(t, df) { // Approximate t-distribution CDF using relationship to F-distribution const x = df / (df + t * t); const p = 0.5 * incompleteBeta(x, df / 2, 0.5); return t >= 0 ? 1 - p : p; } const tValues = model7.beta.map((b, i) => b / stdErrors[i]); const pValues = tValues.map( (t) => 2 * (1 - tCDF(Math.abs(t), df_residual_summary)) ); // R-squared and adjusted R-squared const RSS = RSSResidual; const Rsquared = 1 - RSS / TSS; const adjRsquared = 1 - RSS / df_residual_summary / (TSS / (n - 1)); // Residual standard error const residualStdError = Math.sqrt(sigma2); // F-statistic for overall model const MSS = TSS - RSS; const Fstatistic = MSS / dfModel / (RSS / df_residual_summary); const pValueF = Fstatistic > 100 ? 0 : 1 - fCDF(Fstatistic, dfModel, df_residual_summary); // Format summary output const summaryLines = []; summaryLines.push("\nCall:"); summaryLines.push( "lm(formula = rating ~ lisaProportion * age * gender, data = df)\n" ); summaryLines.push("Residuals:"); summaryLines.push(` Min 1Q Median 3Q Max `); summaryLines.push( `${resMin.toFixed(4).padStart(7)} ${res1Q .toFixed(4) .padStart(7)} ${resMedian.toFixed(4).padStart(7)} ${res3Q .toFixed(4) .padStart(7)} ${resMax.toFixed(4).padStart(7)} \n` ); summaryLines.push("Coefficients:"); summaryLines.push( " Estimate Std. Error t value" ); for (let i = 0; i < nCoef; i++) { const name = coefNames[i].padEnd(45); const est = model7.beta[i].toFixed(6).padStart(9); const se = stdErrors[i].toFixed(6).padStart(11); const tv = tValues[i].toFixed(3).padStart(8); summaryLines.push(`${name}${est}${se}${tv}`); } summaryLines.push(" Pr(>|t|) "); for (let i = 0; i < nCoef; i++) { const name = coefNames[i].padEnd(45); const pv = formatPValue(pValues[i]).padStart(11); const sig = formatSigCode(pValues[i]); summaryLines.push(`${name}${pv}${sig}`); } summaryLines.push("---"); summaryLines.push( "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n" ); summaryLines.push( `Residual standard error: ${residualStdError.toFixed( 1 )} on ${df_residual_summary} degrees of freedom` ); summaryLines.push( `Multiple R-squared: ${Rsquared.toFixed( 5 )},\tAdjusted R-squared: ${adjRsquared.toFixed(5)} ` ); summaryLines.push( `F-statistic: ${Fstatistic.toFixed( 1 )} on ${dfModel} and ${df_residual_summary} DF, p-value: ${formatPValue( pValueF )}\n` ); const df_lisaProp = 1; const df_age = ageUnique.length - 1; const df_gender = genderUnique.length - 1; const df_lisaAge = df_age; const df_lisaGender = df_gender; const df_ageGender = df_age * df_gender; const df_lisaAgeGender = df_age * df_gender; const df_residual = n - (1 + df_lisaProp + df_age + df_gender + df_lisaAge + df_lisaGender + df_ageGender + df_lisaAgeGender); const MS_lisaProp = SS_lisaProp / df_lisaProp; const MS_age = SS_age / df_age; const MS_gender = SS_gender / df_gender; const MS_lisaAge = SS_lisaAge / df_lisaAge; const MS_lisaGender = SS_lisaGender / df_lisaGender; const MS_ageGender = SS_ageGender / df_ageGender; const MS_lisaAgeGender = SS_lisaAgeGender / df_lisaAgeGender; const MS_residual = RSSResidual / df_residual; const F_lisaProp = MS_lisaProp / MS_residual; const F_age = MS_age / MS_residual; const F_gender = MS_gender / MS_residual; const F_lisaAge = MS_lisaAge / MS_residual; const F_lisaGender = MS_lisaGender / MS_residual; const F_ageGender = MS_ageGender / MS_residual; const F_lisaAgeGender = MS_lisaAgeGender / MS_residual; const p_lisaProp = F_lisaProp > 100 ? 0 : 1 - fCDF(F_lisaProp, df_lisaProp, df_residual); const p_age = F_age > 100 ? 0 : 1 - fCDF(F_age, df_age, df_residual); const p_gender = F_gender > 100 ? 0 : 1 - fCDF(F_gender, df_gender, df_residual); const p_lisaAge = 1 - fCDF(F_lisaAge, df_lisaAge, df_residual); const p_lisaGender = 1 - fCDF(F_lisaGender, df_lisaGender, df_residual); const p_ageGender = F_ageGender > 100 ? 0 : 1 - fCDF(F_ageGender, df_ageGender, df_residual); const p_lisaAgeGender = 1 - fCDF(F_lisaAgeGender, df_lisaAgeGender, df_residual); const lines = []; lines.push("Analysis of Variance Table\n"); lines.push("Response: rating"); lines.push( " Df Sum Sq Mean Sq F value Pr(>F) " ); const rows = [ [ "lisaProportion", df_lisaProp, SS_lisaProp, MS_lisaProp, F_lisaProp, p_lisaProp ], ["age", df_age, SS_age, MS_age, F_age, p_age], ["gender", df_gender, SS_gender, MS_gender, F_gender, p_gender], [ "lisaProportion:age", df_lisaAge, SS_lisaAge, MS_lisaAge, F_lisaAge, p_lisaAge ], [ "lisaProportion:gender", df_lisaGender, SS_lisaGender, MS_lisaGender, F_lisaGender, p_lisaGender ], [ "age:gender", df_ageGender, SS_ageGender, MS_ageGender, F_ageGender, p_ageGender ], [ "lisaProportion:age:gender", df_lisaAgeGender, SS_lisaAgeGender, MS_lisaAgeGender, F_lisaAgeGender, p_lisaAgeGender ], ["Residuals", df_residual, RSSResidual, MS_residual, null, null] ]; for (const [name, df, ss, ms, f, p] of rows) { const dfStr = df.toString(); const namePlusSpace = 31 - dfStr.length; // Ensure name+space+df = 32 const nameStr = name.padEnd(namePlusSpace); const ssStr = Math.round(ss).toString().padStart(8); const msStr = ms.toFixed(1).padStart(8); if (f === null) { // Residuals line needs trailing spaces to match F value, Pr(>F), and sig code columns const trailing = " ".repeat(23); // 9 (F value) + 10 (Pr(>F)) + 4 (sig code) lines.push(`${nameStr}${dfStr}${ssStr}${msStr}${trailing}`); } else { const fStr = f.toFixed(3).padStart(9); const pStr = formatPValue(p).padStart(10); const sigCode = formatSigCode(p); lines.push(`${nameStr}${dfStr}${ssStr}${msStr}${fStr}${pStr}${sigCode}`); } } lines.push("---"); lines.push("Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"); return summaryLines.join("\n") + lines.join("\n"); }36:[{"id":"cb8ba3c43c3599e1044dd3a8406212fb23902b7df90b2087da6622b2372c3b96a4cc21615317e033300a1f3848dd3b415109f421274a50b0d3adbb729b07ec24","href":"https://static.observableusercontent.com/files/cb8ba3c43c3599e1044dd3a8406212fb23902b7df90b2087da6622b2372c3b96a4cc21615317e033300a1f3848dd3b415109f421274a50b0d3adbb729b07ec24","name":"all_ratingsbydemo_20200824_115720.tsv","type":"text/tab-separated-values","lastModified":1742773156601.331,"size":8638377},{"id":"c6dd6547ff6dbed0f1c4696520396234241adb8b3c4c7c51bd40480e21acf7f7e8141118175694d6a89fa5d435f71e97e16e3b255803aeedd64e7f49ee150cb5","href":"https://static.observableusercontent.com/files/c6dd6547ff6dbed0f1c4696520396234241adb8b3c4c7c51bd40480e21acf7f7e8141118175694d6a89fa5d435f71e97e16e3b255803aeedd64e7f49ee150cb5","name":"image.png","type":"image/png","lastModified":1742777242504.39,"size":358178}] 30:{"title":"Do male IMDb users dislike Simpson's episodes with more Lisa screentime?","cells":[{"id":0,"value":"# Do male IMDb users dislike Simpson's episodes with more Lisa screentime?\n\n## and is it specifically _young_ men? \n[Max Bo](https://maxbo.me), [Jamie Rosenberg](https://www.linkedin.com/in/jamescrosenberg)\n","pinned":false,"hidden":false,"mode":"md"},{"id":1062,"value":"This notebook investigates 2 hypotheses that my good friend [Jamie Rosenberg](https://www.linkedin.com/in/jamescrosenberg) pitched to me one Friday evening. Fortunately, he had legitimately acquired much Simpsons ratings data from IMDb prior to the removal of their per age and gender cohort ratings on the 2020-08-22.","pinned":false,"hidden":false,"mode":"md"},{"id":1421,"value":"## The Golden Age","pinned":false,"hidden":false,"mode":"md"},{"id":1070,"value":"Firstly, we must ask: _which_ episodes? The Simpsons has a [well documented](https://en.wikipedia.org/wiki/Reception_of_The_Simpsons) \"golden age\", commonly regarded as ending on [episode 203 (S09E25) Natural Born Kissers](https://en.wikipedia.org/wiki/Natural_Born_Kissers):","pinned":false,"hidden":false,"mode":"md"},{"id":940,"value":"viewof goldenAgeFinalEpisodeSequenceNumber = Inputs.range([1, 250], {\n label: \"Final golden age episode sequence number\",\n value: productionCodeToEpisodeSequenceNumber[\"S09E25\"], //‎Natural Born Kissers\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1199,"value":"| field | value |\n| ------ | -------- |\n| Production code | `${cohortAllCohortAverages[goldenAgeFinalEpisodeIndex].productionCode}` |\n| Episode sequence number | `${goldenAgeFinalEpisodeSequenceNumber}` |\n| Title | `${cohortAllCohortAverages[goldenAgeFinalEpisodeIndex].title}` |","pinned":false,"hidden":false,"mode":"md"},{"id":897,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1415,"value":"We can see a fairly obvious partition in quality:","pinned":false,"hidden":false,"mode":"md"},{"id":1021,"value":"split = Plot.plot({\n height: 640,\n marginLeft: 100,\n x: {\n domain: [0, 10],\n label: \"Rating\"\n },\n color: {\n label: \"Rating\",\n type: \"linear\",\n scheme: \"RdYlGn\",\n legend: true,\n domain: [0, 10],\n unknown: \"#ddd\"\n },\n grid: true,\n marks: [\n Plot.frame(),\n Plot.dot(\n cohortAllCohortAverages,\n Plot.dodgeY({\n padding: -1,\n r: 3,\n tip: true,\n channels: {\n Title: (d) => (d.title ? d.title : \"???\"),\n productionCode: \"productionCode\"\n },\n x: \"cohortAverageRating\",\n fill: \"cohortAverageRating\",\n fy: (d, i) =>\n i <= goldenAgeFinalEpisodeIndex ? \"golden age\" : \"post golden age\"\n })\n )\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":51,"value":"## When does Lisa speak?","pinned":false,"hidden":false,"mode":"md"},{"id":1429,"value":"For every Golden Age episode, let's take a look at the script and identify every single instance Lisa speaks:","pinned":false,"hidden":false,"mode":"md"},{"id":40,"value":"Plot.plot({\n y: { grid: true },\n marginLeft: 50,\n height: 640,\n y: {\n ticks: 100\n },\n marks: [\n Plot.frame(),\n Plot.tickX(goldenAgeLines, {\n x: \"timestamp_in_ms\",\n y: \"episode_sequence_number\",\n stroke: \"character\",\n strokeWidth: 3,\n channels: {\n word_count: \"word_count\"\n },\n strokeOpacity: (d) => (d.character === \"Lisa Simpson\" ? 1 : 0.05),\n tip: true\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1554,"value":"And now as a proportion of all words spoken per episode. We see that episode ${exampleDot.episodeSequenceNumber} (${exampleDot.productionCode}), ${exampleDot.title}, is the episode with the most Lisa lines in it (${(exampleDot.lisaProportion * 100).toFixed(1)}%):","pinned":false,"hidden":false,"mode":"md"},{"id":1475,"value":"Plot.plot({\n nice: true,\n height: 640,\n x: {\n grid: true,\n label: \"Proportion of words spoken (%)\",\n percent: true\n },\n y: {\n domain: sortedByLisaProportionLisaEpisodes\n },\n marks: [\n Plot.frame(),\n Plot.barX(\n sortedByLisaGoldenAgeLines,\n Plot.stackX(\n {\n order: null,\n offset: \"expand\"\n },\n Plot.groupY(\n {\n x: \"sum\",\n fillOpacity: (d) => (d[0].character === \"Lisa Simpson\" ? 1 : 0.1)\n },\n {\n fill: \"character\",\n tip: true,\n y: \"episode_sequence_number\",\n x: \"word_count\"\n }\n )\n )\n )\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2256,"value":"From that, let's plot the relationship between ratings and the proportion of words spoken by Lisa. Each dot's y-axis is the the average vote of an entire cohort.","pinned":false,"hidden":false,"mode":"md"},{"id":2275,"value":"For example, episode ${exampleDot.productionCode} (_${exampleDot.title}_), for cohort `${exampleDot.age}` / `${exampleDot.gender}` has an average rating of ⬤ ${exampleDot.cohortAverageRating.toFixed(1)}. The regression line is NOT against these average cohort ratings, but against the \"synthetic\" data points that would otherwise be at the specified rating levels. Each dot is scaled by total vote count for that episode in that cohort.","pinned":false,"hidden":false,"mode":"md"},{"id":674,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2244,"value":"## Stats","pinned":false,"hidden":false,"mode":"md"},{"id":2246,"value":"`model <- lm(rating ~ lisaProportion * age * gender, data = df)`","pinned":false,"hidden":false,"mode":"md"},{"id":1163,"value":"anova(\n exploded.filter(\n (d) => d.lisaProportion && d.age !== \"all\" && d.gender !== \"all\"\n )\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2483,"value":"Male IMDb users **do not** appear to dislike Lisa episodes more specifically because of age. _However_, there is evidence of a gender difference in the relationship between Lisa screentime and ratings: male cohorts rate higher-Lisa episodes lower relative to female cohorts. The model does not support the stronger claim that young men uniquely drive this effect.","pinned":false,"hidden":false,"mode":"md"},{"id":1349,"value":"---\n\n## Implementation","pinned":false,"hidden":false,"mode":"md"},{"id":1993,"value":"### Golden age","pinned":false,"hidden":false,"mode":"md"},{"id":1940,"value":"Style the episode slider (well, all sliders)","pinned":false,"hidden":false,"mode":"md"},{"id":2081,"value":"goldenAgeFinalEpisodeIndex = goldenAgeFinalEpisodeSequenceNumber - 1","pinned":true,"hidden":false,"mode":"ojs"},{"id":994,"value":"","pinned":true,"hidden":false,"mode":"md"},{"id":911,"value":"import { colorContrast } from \"@observablehq/plot-colorcontrast-custom-transform\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":998,"value":"import { settle } from \"@mjbo/settling-input\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1369,"value":"When did the user specify the \"golden age\" cutoff?","pinned":false,"hidden":false,"mode":"md"},{"id":1003,"value":"settledgoldenAgeFinalEpisodeSequenceNumber = settle(viewof goldenAgeFinalEpisodeSequenceNumber)","pinned":true,"hidden":false,"mode":"ojs"},{"id":972,"value":"settledGoldenAgeFinalEpisodeIndex = settledgoldenAgeFinalEpisodeSequenceNumber - 1","pinned":true,"hidden":false,"mode":"ojs"},{"id":1995,"value":"### Lisa proportion","pinned":false,"hidden":false,"mode":"md"},{"id":1373,"value":"Pull in all the lines:","pinned":false,"hidden":false,"mode":"md"},{"id":4,"value":"import { metaLines } from \"@enjalot/simpsons-data\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1970,"value":"Drop and rename fields:","pinned":false,"hidden":false,"mode":"md"},{"id":1427,"value":"lines = metaLines.map((d) => {\n return {\n episode_sequence_number: d.episode_id,\n timestamp_in_ms: d.timestamp_in_ms,\n speaking_line: d.speaking_line,\n character: d.raw_character_text,\n word_count: d.word_count\n };\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":1384,"value":"Script episode sequence numbers are 1-indexed:","pinned":false,"hidden":false,"mode":"md"},{"id":1375,"value":"[...new Set(lines.map((d) => d.episode_sequence_number))].toSorted(\n (a, b) => a - b\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1386,"value":"Filter down to just golden age lines: ","pinned":false,"hidden":false,"mode":"md"},{"id":74,"value":"goldenAgeLines = lines.filter(\n (d) =>\n d.episode_sequence_number <= settledgoldenAgeFinalEpisodeSequenceNumber &&\n d.speaking_line &&\n !isNaN(d.word_count) &&\n d.word_count < 100_000 // we have 3 weird outliers with massive wordcounts, see below:\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1362,"value":"lines.filter((d) => d.word_count > 100_000)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1539,"value":"This is a bit of a hack, but this is so that Lisa always appears first in [Stack transform ordering](https://observablehq.com/plot/transforms/stack).","pinned":false,"hidden":false,"mode":"md"},{"id":1536,"value":"sortedByLisaGoldenAgeLines = goldenAgeLines.sort((a, b) => {\n const aIsLisa = a.character === \"Lisa Simpson\";\n const bIsLisa = b.character === \"Lisa Simpson\";\n return bIsLisa - aIsLisa;\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":1948,"value":"For each episode number, compute the proportion that Lisa speaks out of all the words.","pinned":false,"hidden":false,"mode":"md"},{"id":599,"value":"episodeSequenceNumberToLisaProportion = {\n const epToLines = {};\n const epToLisaLines = {};\n const epToProportion = {};\n\n for (const line of goldenAgeLines) {\n epToLines[line.episode_sequence_number] =\n (epToLines[line.episode_sequence_number] ?? 0) + line.word_count;\n\n if (line.character === \"Lisa Simpson\") {\n epToLisaLines[line.episode_sequence_number] =\n (epToLisaLines[line.episode_sequence_number] ?? 0) + line.word_count;\n }\n }\n\n for (const ep of Object.keys(epToLines)) {\n epToProportion[ep] = epToLisaLines[ep] / epToLines[ep];\n }\n\n return epToProportion;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1991,"value":"We need this to sort the y-domain of all our charts:","pinned":false,"hidden":false,"mode":"md"},{"id":1466,"value":"sortedByLisaProportionLisaEpisodes = Object.entries(\n episodeSequenceNumberToLisaProportion\n)\n .toSorted((a, b) => b[1] - a[1])\n .map((d) => Number(d[0]))","pinned":true,"hidden":false,"mode":"ojs"},{"id":636,"value":"### IMDb data acquired through legitimate and legal means","pinned":false,"hidden":false,"mode":"md"},{"id":2425,"value":"ratings = FileAttachment(\"all_ratingsbydemo_20200824_115720.tsv\").tsv({\n typed: true\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":647,"value":"Our IMDb data set doesn't have episode numbers, only production codes. So we need to derive a production code -> episode number mapping:","pinned":false,"hidden":false,"mode":"md"},{"id":649,"value":"productionCodeToEpisodeSequenceNumber = {\n const mapping = {};\n [...new Set(ratings.map((d) => d.seNum))].toSorted().forEach((seNum, i) => {\n mapping[seNum] = i + 1;\n });\n\n return mapping;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2013,"value":"We also construct a lookup table of production code to episode metadata (we're mainly interested in the episode title):","pinned":false,"hidden":false,"mode":"md"},{"id":925,"value":"import { simpsons as episodeMetadata } from \"@observablehq/plot-simpsons-ratings\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2319,"value":"episodeSequenceNumberToMetadata = {\n const acc = {};\n\n for (const d of episodeMetadata) {\n acc[d.id] = d;\n }\n\n return acc;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2061,"value":"### Join","pinned":false,"hidden":false,"mode":"md"},{"id":2016,"value":"Now we join our IMDb ratings, metadata and the information we gathered from the scripts. We also:\n\n1. Derive a \"cohort\" which we use to key a bunch of computations later on\n1. Filter the age cohorts to those that we care about","pinned":false,"hidden":false,"mode":"md"},{"id":2174,"value":"AGE_DOMAIN = new Set([\"all\", \"aged_under_18\", \"aged_18_29\", \"aged_30_44\", \"aged_45_plus\"])","pinned":true,"hidden":true,"mode":"ojs"},{"id":2177,"value":"GENDER_DOMAIN = new Set([\"males\", \"all\", \"females\"])","pinned":true,"hidden":true,"mode":"ojs"},{"id":668,"value":"tidy = ratings\n .map((d) => {\n return {\n ...d,\n productionCode: d.seNum,\n episodeSequenceNumber: productionCodeToEpisodeSequenceNumber[d.seNum]\n };\n })\n .map((d) => ({\n episodeSequenceNumber: d.episodeSequenceNumber,\n productionCode: d.productionCode,\n ...parseSeasonAndEpisode(d.productionCode),\n lisaProportion:\n episodeSequenceNumberToLisaProportion[\n productionCodeToEpisodeSequenceNumber[d.productionCode]\n ],\n rating: d.rating,\n votes: d.votes,\n title: episodeSequenceNumberToMetadata[d.episodeSequenceNumber]?.title,\n gender: d.cat_b === \"gender\" ? d.cat_c : d.cat_b,\n age: d.cat_b === \"gender\" ? \"all\" : d.cat_c,\n cohort:\n d.episodeSequenceNumber + \"|\" + d.cat_a + \"|\" + d.cat_b + \"|\" + d.cat_c\n }))\n .filter((d) => AGE_DOMAIN.has(d.age) && GENDER_DOMAIN.has(d.gender))","pinned":true,"hidden":false,"mode":"ojs"},{"id":686,"value":"cohortToAverageRatingAndTotalVotesPerEp = {\n const cohortTotalVotes = {};\n const cohortTotalScore = {};\n const cohortAverageRating = {};\n\n for (const d of tidy) {\n cohortTotalVotes[d.cohort] = (cohortTotalVotes[d.cohort] ?? 0) + d.votes;\n cohortTotalScore[d.cohort] =\n (cohortTotalScore[d.cohort] ?? 0) + d.votes * d.rating;\n }\n\n for (const cohort of Object.keys(cohortTotalScore)) {\n cohortAverageRating[cohort] =\n cohortTotalScore[cohort] / cohortTotalVotes[cohort];\n }\n\n return { cohortAverageRating, cohortTotalVotes };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":906,"value":"function parseSeasonAndEpisode(productionCode) {\n const match = productionCode.match(/S(\\d+)E(\\d+)/i);\n if (match) {\n return {\n season: parseInt(match[1]),\n episode: parseInt(match[2])\n };\n }\n throw new Error(\"failed parse\");\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2057,"value":"Now we join even more episode in, which is the average rating for that episode for that specific cohort:","pinned":false,"hidden":false,"mode":"md"},{"id":2260,"value":"exampleDot = cohortAverages\n .filter((d) => d.age === \"aged_18_29\" && d.gender === \"males\")\n .toSorted((a, b) => b.lisaProportion - a.lisaProportion)[0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":706,"value":"cohortAverages = {\n const intermediate = tidy.map((d) => ({\n ...d,\n cohortAverageRating:\n cohortToAverageRatingAndTotalVotesPerEp.cohortAverageRating[d.cohort],\n cohortTotalVotes:\n cohortToAverageRatingAndTotalVotesPerEp.cohortTotalVotes[d.cohort]\n }));\n\n return _.uniqBy(intermediate, \"cohort\");\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2111,"value":"`exploded` creates \"synthetic\" datapoints based on the number of votes a certain cohort gave for a certain rating. We use this to generate our linear regression lines.","pinned":false,"hidden":false,"mode":"md"},{"id":824,"value":"exploded = {\n const acc = [];\n\n for (const d of tidy) {\n for (let i = 0; i < d.votes; ++i) {\n acc.push(d);\n }\n }\n\n return acc;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":902,"value":"cohortAllCohortAverages = cohortAverages\n .filter((d) => d.gender === \"all\" && d.age === \"all\")\n .filter((d) => d.cohortAverageRating)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2222,"value":"### Statistical analysis","pinned":false,"hidden":false,"mode":"md"},{"id":2212,"value":"I built this linear model in R and then asked Claude to, character for character, replicate the output with equivalent JavaScript implementation.","pinned":false,"hidden":false,"mode":"md"},{"id":2214,"value":"$34","pinned":false,"hidden":false,"mode":"md"},{"id":2206,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}