MentDB WeakMQL Recipes
← innov-ai.com All profiles Portal
Profiles / Engineering & data
SB

Statistician & biostatistician

Confidence intervals, sample size, ANOVA, chi-square and non-parametric tests.

10 recipes on this page · open in the interactive view

Recipes
  1. Confidence interval and sample size
  2. ANOVA and Kruskal-Wallis on three diets
  3. Clinical trial 2x2 table
  4. Bayesian A/B test with Beta posteriors
  5. Before / after treatment: paired t-test or Wilcoxon
  6. z-test with known sigma, power and sample size
  7. Summary of a 5-point Likert survey item
  8. Which mean? Arithmetic, geometric and harmonic in practice
  9. The Student t density from the Gamma and Beta functions
  10. Mincer wage equation from a survey table in SQL
executed on MentDB

Confidence interval and sample size

A pilot study measures the drop of systolic pressure on 10 patients. The script gives the 95 % confidence interval of the mean and the number of patients needed for a margin of 1.5 mmHg.

Pipeline
01 · intervalStudent confidence intervalstat conf_interval_mean→02 · quantilez for 95 %stat dist_inv→03 · sizen = (z s / E)^2 rounded upmath ceil
Script
admin
#1. Pilot study: drop of systolic pressure in mmHg after 8 weeks, 10 patients;
-> "[d]" "[8, 12, 5, 15, 10, 7, 13, 9, 11, 6]";
json load "rx_st_ci" (stat conf_interval_mean [d] 0.95);
#2. Sample size for a margin of 1.5 mmHg at 95 percent: n = z s / E squared;
-> "[z]" (stat dist_inv "normal" "[0, 1]" 0.975);
-> "[n]" (math ceil (math calc (concat "([z]*" (stat stddev [d]) "/1.5)^2")));
concat "mean = " (json select "rx_st_ci" /mean) " mmHg, 95% CI [" (math round (json select "rx_st_ci" /lower) 2) ", " (math round (json select "rx_st_ci" /upper) 2) "], margin " (math round (json select "rx_st_ci" /margin) 2) " -> " [n] " patients needed for a margin of 1.5 mmHg"
Result
mentdb
mean = 9.6 mmHg, 95% CI [7.31, 11.89], margin 2.29 -> 18 patients needed for a margin of 1.5 mmHg
Go further
  • Use the Student quantile with stat dist_inv "student" and iterate on n
  • Add the power of a test: n = ((z_alpha + z_beta) s / delta)^2
  • Check the hypothesis of the paired difference with stat ttest_1
executed on MentDB

ANOVA and Kruskal-Wallis on three diets

Weight loss is measured on three diets of six subjects. The script runs a one-way ANOVA with its effect size and the non-parametric Kruskal-Wallis test as a check.

Pipeline
01 · anovaOne-way ANOVAstat anova→02 · ranksKruskal-Wallis teststat kruskal→03 · reportF, p, eta2, Hmath round_sig
Script
admin
#1. Weight loss in kg with three diets;
-> "[g]" "[[3.1, 4.2, 2.8, 3.9, 3.5, 4.0], [5.2, 4.8, 6.1, 5.5, 4.9, 5.8], [3.8, 4.5, 4.1, 3.6, 4.9, 4.4]]";
#2. Parametric one-way ANOVA;
json load "rx_st_a" (stat anova [g]);
#3. Non-parametric check: Kruskal-Wallis on the ranks;
json load "rx_st_k" (stat kruskal [g]);
concat "ANOVA F(" (json select "rx_st_a" /df_between) ", " (json select "rx_st_a" /df_within) ") = " (math round (json select "rx_st_a" /f) 2) ", p = " (math round_sig (json select "rx_st_a" /p_value) 3) ", eta2 = " (math round (json select "rx_st_a" /eta2) 3) " ; Kruskal-Wallis H = " (math round (json select "rx_st_k" /h) 2) ", p = " (math round_sig (json select "rx_st_k" /p_value) 3)
Result
mentdb
ANOVA F(2, 15) = 18.92, p = 0.0000791, eta2 = 0.716 ; Kruskal-Wallis H = 12.07, p = 0.00240
Go further
executed on MentDB

Clinical trial 2x2 table

A trial counts the events in a treatment group and a placebo group of 100 patients each. The script runs the chi-square and Fisher exact tests and computes the relative risk, the absolute risk reduction and the number needed to treat.

Pipeline
01 · chi2Chi-square of independencestat chi2_independence→02 · exactFisher exact teststat fisher_exact→03 · effectRR, ARR, NNTmath calc
Script
admin
#1. Trial: events / no event, treatment 12 / 88, placebo 27 / 73;
-> "[t]" "[[12, 88], [27, 73]]";
#2. Chi-square test of independence and Fisher exact test;
json load "rx_st_c" (stat chi2_independence [t]);
json load "rx_st_f" (stat fisher_exact [t]);
#3. Relative risk, absolute risk reduction and number needed to treat;
-> "[rr]" (math calc "(12/100)/(27/100)");
-> "[arr]" (math calc "27/100-12/100");
concat "chi2 = " (math round (json select "rx_st_c" /chi2) 2) ", p = " (math round_sig (json select "rx_st_c" /p_value) 3) " ; Fisher p = " (math round_sig (json select "rx_st_f" /p_value) 3) ", OR = " (math round (json select "rx_st_f" /odds_ratio) 3) " ; RR = " (math round [rr] 3) ", ARR = " (math round (math calc "100*[arr]") 1) " %, NNT = " (math ceil (math calc "1/[arr]"))
Result
mentdb
chi2 = 7.17, p = 0.00743 ; Fisher p = 0.0118, OR = 0.369 ; RR = 0.444, ARR = 15 %, NNT = 7
Go further
executed on MentDB

Bayesian A/B test with Beta posteriors

Conversions 120/1000 against 145/1000 with uniform priors give Beta posteriors. The script computes the probability that B is better than A by a double integral of the joint posterior (normalised with lgamma) and checks it in one dimension with the incomplete beta function.

Pipeline
01 · normln B(a, b) with lgammalgamma(x)→02 · 2dDouble integral over y > xnum integrate2→03 · 1dCheck with I_x(a, b)beta_inc(a, b, x)
Script
admin
#1. A/B test: A 120 conversions out of 1000, B 145 out of 1000, uniform priors: Beta posteriors;
-> "[a1]" 121; -> "[b1]" 881; -> "[a2]" 146; -> "[b2]" 856;
#2. Log normalising constants with lgamma: ln B of a,b = lgamma a + lgamma b - lgamma a+b;
-> "[lB1]" (math calc "lgamma([a1])+lgamma([b1])-lgamma([a1]+[b1])");
-> "[lB2]" (math calc "lgamma([a2])+lgamma([b2])-lgamma([a2]+[b2])");
#3. P of pB > pA: double integral of the joint posterior over y > x;
-> "[P]" (num integrate2 (concat "exp(120*ln(x)+880*ln(1-x)-(" [lB1] ")+145*ln(y)+855*ln(1-y)-(" [lB2] "))") "x" 0.05 0.25 "y" "x" 0.3);
#4. Check in one dimension with the incomplete beta function;
-> "[P1]" (num integrate (concat "exp(120*ln(x)+880*ln(1-x)-(" [lB1] "))*(1-beta_inc(146, 856, x))") "x" 0.05 0.25);
concat "P(B better than A) = " (math round [P] 4) " (1D check " (math round [P1] 4) "), posterior means " (math round (math calc "[a1]/([a1]+[b1])") 4) " and " (math round (math calc "[a2]/([a2]+[b2])") 4)
Result
mentdb
P(B better than A) = 0.9502 (1D check 0.9502), posterior means 0.1208 and 0.1457
Go further
executed on MentDB

Before / after treatment: paired t-test or Wilcoxon

The systolic pressure of 12 patients is measured before and after 8 weeks of treatment. Test the normality of the differences (Shapiro-Wilk), run the paired t-test and the Wilcoxon signed-rank test, and give the effect size Cohen dz. Here the differences are not normal (p < 0.05), so the Wilcoxon test is the one to report.

Pipeline
01 · differencesBefore - aftervector sub→02 · normalityShapiro-Wilk on dstat shapiro→03 · parametricPaired t-teststat ttest_paired→04 · rankWilcoxon signed-rankstat wilcoxon→05 · effectCohen dzstat stddev
Script
admin
#1. Systolic blood pressure in mmHg of 12 patients before and after 8 weeks of treatment;
-> "[before]" "[152, 148, 160, 145, 158, 151, 163, 149, 155, 147, 159, 150]";
-> "[after]" "[141, 145, 148, 140, 150, 147, 151, 146, 144, 145, 149, 148]";
#2. Differences and their normality, which the paired t-test assumes;
-> "[d]" (vector sub [before] [after]);
json load "rx_sw" (stat shapiro [d]);
#3. Paired t-test and its non-parametric counterpart;
json load "rx_tp" (stat ttest_paired [before] [after]);
json load "rx_wx" (stat wilcoxon [before] [after]);
#4. Effect size: Cohen dz = mean of d / sd of d;
-> "[dz]" (/ (stat mean [d]) (stat stddev [d]));
concat "mean drop " (math round (stat mean [d]) 2) " mmHg | Shapiro W = " (math round (json select "rx_sw" "/w") 3) ", p = " (math round (json select "rx_sw" "/p_value") 3) " | paired t = " (math round (json select "rx_tp" "/t") 3) ", df = " (json select "rx_tp" "/df") ", p = " (math round_sig (json select "rx_tp" "/p_value") 3) " | Wilcoxon W = " (json select "rx_wx" "/w") ", p = " (math round_sig (json select "rx_wx" "/p_value") 3) " (" (json select "rx_wx" "/method") ") | Cohen dz = " (math round [dz] 2)
Result
mentdb
mean drop 6.92 mmHg | Shapiro W = 0.848, p = 0.035 | paired t = 5.812, df = 11, p = 0.000117 | Wilcoxon W = 0, p = 0.000488 (permutation) | Cohen dz = 1.68
Go further
executed on MentDB

z-test with known sigma, power and sample size

A standardised reading test has a known sigma of 15 and a national mean of 100. A school tested 12 pupils. The statistician runs the z-test, computes the power to detect a 5-point shift with this n, and the n needed for 80 % power, all from the error function.

Pipeline
01 · testz and p-valuestat ztest→02 · critical valuez = sqrt2 erf_inv(1 - alpha)math erf_inv→03 · powerPhi through erfcmath erfc→04 · sample sizen = ((z_a + z_b) sigma / delta)^2math ceil
Script
admin
#1. Scores of the class, H0: mu = 100, sigma = 15;
-> "[x]" "[112, 98, 105, 121, 93, 108, 117, 101, 99, 110, 125, 104]";
json load "rx_p4g_z" (stat ztest [x] 100 15);
-> "[t]" (concat "mean = " (math round (stat mean [x]) 2) ", z = " (math round (json select "rx_p4g_z" /z) 3) ", p = " (math round (json select "rx_p4g_z" /p_value) 4));
json unload "rx_p4g_z";
#2. Critical values: Phi^-1 of p = sqrt2 . erf_inv of 2p - 1;
-> "[za]" (* (math sqrt 2) (math erf_inv 0.95));
-> "[zb]" (* (math sqrt 2) (math erf_inv 0.6));
#3. Power for delta = 5 with n = 12, Phi of x = 0.5 . erfc of -x / sqrt2;
-> "[k]" (- (/ (* 5 (math sqrt 12)) 15) [za]);
-> "[pw]" (* 0.5 (math erfc (/ (- 0 [k]) (math sqrt 2))));
#4. n for 80 % power;
-> "[n]" (math ceil (math power (/ (* (+ [za] [zb]) 15) 5) 2));
concat [t] [_n_] "z 0.975 = " (math round [za] 4) ", z 0.80 = " (math round [zb] 4) [_n_] "power for +5 points with n = 12: " (math round (* 100 [pw]) 1) " %" [_n_] "n for 80 % power: " [n]
Result
mentdb
mean = 107.75, z = 1.79, p = 0.0735
z 0.975 = 1.96, z 0.80 = 0.8416
power for +5 points with n = 12: 21 %
n for 80 % power: 71
Go further
  • Compare with stat ttest_1 when sigma is estimated from the sample.
  • Plot the power curve against n with plot chart.
executed on MentDB

Summary of a 5-point Likert survey item

An item of a satisfaction survey is ordinal (1 = strongly disagree to 5 = strongly agree). The statistician reports what fits ordinal data: the frequency table, the mode, the median, the top-2-box and the bottom-2-box, with the mean only for information.

Pipeline
01 · tableCounts and cumulative frequenciesstat frequency→02 · modeMost frequent answerstat mode→03 · medianMedian answerstat median→04 · boxesTop-2 and bottom-2 boxjson select→05 · meanMean for informationmath avg
Script
admin
#1. Answers of 24 respondents;
-> "[a]" "[4, 5, 3, 4, 4, 2, 5, 4, 3, 1, 4, 5, 4, 3, 2, 4, 5, 4, 3, 4, 5, 2, 4, 3]";
json load "rx_p4g_lk" (stat frequency [a]);
#2. Frequency per answer;
-> "[tab]" "";
for (-> "[i]" 0) (< [i] (json count "rx_p4g_lk" /)) (++ "[i]") {
  -> "[tab]" (concat [tab] (json select "rx_p4g_lk" (concat "/[" [i] "]/value")) ": " (math round (* 100 (json select "rx_p4g_lk" (concat "/[" [i] "]/frequency"))) 1) " %  ");
};
#3. Bottom-2 box = cumulative frequency of answer 2, top-2 = 1 - cumulative of 3;
-> "[b2]" (json select "rx_p4g_lk" "/[1]/cumulative_frequency");
-> "[t2]" (- 1 (json select "rx_p4g_lk" "/[2]/cumulative_frequency"));
json unload "rx_p4g_lk";
concat [tab] [_n_] "mode = " (stat mode [a]) ", median = " (stat median [a]) ", mean = " (math round (math avg [a]) 2) [_n_] "top-2-box = " (math round (* 100 [t2]) 1) " %, bottom-2-box = " (math round (* 100 [b2]) 1) " %, net score = " (math round (* 100 (- [t2] [b2])) 1)
Result
mentdb
1: 4.2 %  2: 12.5 %  3: 20.8 %  4: 41.7 %  5: 20.8 %
mode = [4], median = 4, mean = 3.63
top-2-box = 62.5 %, bottom-2-box = 16.7 %, net score = 45.8
Go further
executed on MentDB

Which mean? Arithmetic, geometric and harmonic in practice

Three classic cases where the arithmetic mean misleads: the average yearly growth of an investment (geometric mean), the P/E ratio of an equal-amount portfolio (harmonic mean), and the population versus sample variance of a full census (divide by n).

Pipeline
01 · arithmeticSimple averagemath avg→02 · geometricMean growth factorstat geo_mean→03 · harmonicPortfolio P/Estat harm_mean→04 · variancePopulation vs samplestat pvariance
Script
admin
#1. Yearly growth factors: +25 %, -20 %, +15 %, +5 %;
-> "[g]" "[1.25, 0.80, 1.15, 1.05]";
-> "[ga]" (math avg [g]);
-> "[gg]" (stat geo_mean [g]);
#2. Same amount invested in 4 shares with these P/E ratios;
-> "[pe]" "[12, 18, 30, 45]";
#3. Census of the 6 branches of a company: staff per branch;
-> "[s]" "[14, 22, 9, 31, 18, 26]";
concat "growth: arithmetic " (math round (* 100 (- [ga] 1)) 2) " %/yr is wrong, geometric " (math round (* 100 (- [gg] 1)) 2) " %/yr, check final value " (math round (math power [gg] 4) 4) " = " (math round (* 1.25 (* 0.8 (* 1.15 1.05))) 4) [_n_] "portfolio P/E: harmonic " (math round (stat harm_mean [pe]) 2) " (arithmetic would say " (math avg [pe]) ")" [_n_] "census variance /n = " (math round (stat pvariance [s]) 3) ", sample variance /(n-1) = " (math round (stat variance [s]) 3)
Result
mentdb
growth: arithmetic 6.25 %/yr is wrong, geometric 4.83 %/yr, check final value 1.2075 = 1.2075
portfolio P/E: harmonic 20.57 (arithmetic would say 26.25)
census variance /n = 53.667, sample variance /(n-1) = 64.4
Go further
  • Show the inequality harmonic <= geometric <= arithmetic on any positive series.
  • Use stat weighted_mean when the amounts invested differ.
executed on MentDB

The Student t density from the Gamma and Beta functions

To teach where the t density comes from, the script evaluates f(t) = Gamma((v+1)/2) / (sqrt(v pi) Gamma(v/2)) (1 + t^2/v)^(-(v+1)/2) with the Gamma function, then with the Beta form, then with lngamma for v = 200 (where Gamma overflows), and compares with the built-in law.

Pipeline
01 · Gamma formGamma ratiomath gamma→02 · Beta form1 / (sqrt v B(1/2, v/2))math beta→03 · large vexp of lngamma differencesmath lngamma→04 · referenceBuilt-in densitystat dist_pdf
Script
admin
#1. v = 5 degrees of freedom, t = 1.5;
-> "[v]" 5; -> "[t]" 1.5;
-> "[k]" (math power (+ 1 (/ (* [t] [t]) [v])) (- 0 (/ (+ [v] 1) 2)));
-> "[f1]" (* (/ (math gamma (/ (+ [v] 1) 2)) (* (math sqrt (* [v] (math pi))) (math gamma (/ [v] 2)))) [k]);
#2. Same with the Beta function;
-> "[f2]" (/ [k] (* (math sqrt [v]) (math beta 0.5 (/ [v] 2))));
#3. v = 200: Gamma of 100.5 is huge, use lngamma;
-> "[c200]" (math exp (- (math lngamma 100.5) (math lngamma 100)));
-> "[f200]" (* (/ [c200] (math sqrt (* 200 (math pi)))) (math power (+ 1 (/ 2.25 200)) -100.5));
concat "v = 5: Gamma form " (math round [f1] 8) ", Beta form " (math round [f2] 8) ", stat dist_pdf " (math round (stat dist_pdf "student" "[5]" 1.5) 8) [_n_] "v = 200: lngamma form " (math round [f200] 6) ", stat dist_pdf " (math round (stat dist_pdf "student" "[200]" 1.5) 6) ", normal " (math round (stat dist_pdf "normal" "[0, 1]" 1.5) 6)
Result
mentdb
v = 5: Gamma form 0.12451734, Beta form 0.12451734, stat dist_pdf 0.12451734
v = 200: lngamma form 0.129445, stat dist_pdf 0.129445, normal 0.129518
Go further
executed on MentDB

Mincer wage equation from a survey table in SQL

Survey records sit in a database table (embedded H2 here). The statistician loads log hourly wage against years of education and experience straight from a SELECT into a multiple regression, and reports the coefficients, their standard errors and the adjusted R2.

Pipeline
01 · connectionRegister the databasecm set h2_embedded→02 · dataTable of respondentssql dml→03 · modelLoad ln wage ~ educ + exper from SQLpa rm load→04 · estimatesCoefficients and SEpa rm estimate_regression_parameters_standard_errors→05 · fitAdjusted R2pa rm calculate_adjusted_r_squared
Script
admin
#1. Survey database, embedded H2 file under tmp;
cm set "tuto_p4g_svy" {execute "db.h2.embedded.config.get"
	"[database]" "./tmp/tuto_p4g_svy/survey"
	"[user]" "tuto"
	"[password]" "tuto_pwd"
;};
sql connect "tuto_p4g_s2" {cm get "tuto_p4g_svy";};
sql dml "tuto_p4g_s2" "create table resp (educ double, exper double, lnwage double)";
sql dml "tuto_p4g_s2" "insert into resp values (10, 5, 2.45), (12, 3, 2.52), (12, 10, 2.78), (14, 2, 2.70), (16, 4, 2.98), (16, 12, 3.21), (11, 20, 2.80), (18, 6, 3.25), (13, 8, 2.80), (15, 15, 3.12), (12, 25, 2.95), (17, 9, 3.20)";
sql disconnect "tuto_p4g_s2";
#2. Multiple regression straight from the query, unused X fields empty;
pa rm load "rx_p4g_wage" "tuto_p4g_svy" "educ" "exper" "" "" "" "lnwage" "select educ, exper, lnwage from resp";
json load "rx_p4g_b" (pa rm estimate_regression_parameters "rx_p4g_wage");
json load "rx_p4g_se" (pa rm estimate_regression_parameters_standard_errors "rx_p4g_wage");
#3. Return to one year of schooling = exp of b1 - 1;
-> "[b1]" (json select "rx_p4g_b" "/[1]");
-> "[r]" (concat "ln w = " (math round (json select "rx_p4g_b" "/[0]") 3) " + " (math round [b1] 4) " educ + " (math round (json select "rx_p4g_b" "/[2]") 4) " exper" [_n_] "SE: " (math round (json select "rx_p4g_se" "/[0]") 3) ", " (math round (json select "rx_p4g_se" "/[1]") 4) ", " (math round (json select "rx_p4g_se" "/[2]") 4) ", t educ = " (math round (/ [b1] (json select "rx_p4g_se" "/[1]")) 2) [_n_] "adjusted R2 = " (math round (pa rm calculate_adjusted_r_squared "rx_p4g_wage") 4) ", return to schooling = " (math round (* 100 (- (math exp [b1]) 1)) 2) " % per year");
#4. Close and remove the temporary database;
json unload "rx_p4g_b"; json unload "rx_p4g_se";
pa rm close "rx_p4g_wage";
cm remove "tuto_p4g_svy";
file delete "tmp/tuto_p4g_svy";
[r]
Result
mentdb
ln w = 1.296 + 0.1018 educ + 0.0195 exper
SE: 0.087, 0.0057, 0.0021, t educ = 17.93
adjusted R2 = 0.9695, return to schooling = 10.71 % per year
Go further

Also in Engineering & data

TR
Telecom & radar engineer
Link budgets in dB, radar equation, detection thresholds and spectra of modulated signals.
10 recipes
CE
Control engineer
Transfer functions, margins, PID and state feedback: complex numbers, CAS, ODE and matrices.
10 recipes
RE
Robotics engineer
Kinematics, Jacobians and trajectories: matrices, vectors, non-linear systems and splines.
10 recipes
EE
Electrical engineer
Impedances, three-phase power and cable sizing with complex numbers and units.
10 recipes
ES
Electronics & signal engineer
Filter design, spectra, peak detection and SNR on sampled signals.
10 recipes
IS
IoT & sensor engineer
Sensor calibration, smoothing, drift detection and measurement uncertainty.
10 recipes
MC
Mechanical & civil engineer
Deflection, stresses, safety factors and concrete strength statistics.
10 recipes
IQ
Industrial quality engineer
SPC control limits, process capability and tests on machines.
10 recipes
DS
Data scientist
Regression, correlation, feature scaling and drift tests.
10 recipes
DQ
Data quality analyst
Outliers, format checks, normalization and duplicate detection.
10 recipes