{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# IntroStat Week 6\n", "\n", "Welcome to the 6th lecture in IntroStat\n", "\n", "During the lectures we will present both slides and notebooks. \n", "\n", "This is the notebook used in the lecture in week 6.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as stats\n", "import statsmodels.api as sm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Nutrition study 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First visualise the data :)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Enter data\n", "A = np.array([7.53, 7.48, 8.08, 8.09, 10.15, 8.40, 10.88, 6.13, 7.90])\n", "B = np.array([9.21, 11.51, 12.79, 11.85, 9.97, 8.79, 9.69, 9.68, 9.19])\n", "\n", "plt.hist(A, alpha=0.5)\n", "plt.hist(B, alpha=0.5)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets calculate mean and sample standard deviation for each group" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mean_A = A.mean()\n", "mean_B = B.mean()\n", "print([mean_A, mean_B])\n", "\n", "s_A = A.std(ddof=1)\n", "s_B = B.std(ddof=1)\n", "print([s_A, s_B])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the mean values as well" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(A, alpha=0.5)\n", "plt.hist(B, alpha=0.5)\n", "plt.axvline(mean_A, linestyle='-', color=\"blue\", ymin=0, ymax=1)\n", "plt.axvline(mean_B, linestyle='-', color=\"orange\", ymin=0, ymax=1)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great!\n", "But are they the same?\n", "\n", "What else can we do?\n", "\n", "lets calculate standard error (se) for each group mean and add the confidence intervals " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_A = len(A) \n", "n_B = len(B)\n", "se_A = s_A/np.sqrt(n_A)\n", "se_B = s_B/np.sqrt(n_B)\n", "print([se_A, se_B])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Need to choose alpha - we choose alpha = 0.05 \n", "t_0975 = stats.t.ppf(0.975, df=n_A-1) # same for B since n_A = n_B\n", "\n", "CI_A = [mean_A - t_0975 * se_A, mean_A + t_0975 * se_A]\n", "CI_B = [mean_B - t_0975 * se_B, mean_B + t_0975 * se_B]\n", "print(CI_A)\n", "print(CI_B)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(A, alpha=0.5)\n", "plt.hist(B, alpha=0.5)\n", "plt.axvline(mean_A, linestyle='-', color=\"blue\", ymin=0, ymax=1)\n", "plt.axvline(mean_B, linestyle='-', color=\"orange\", ymin=0, ymax=1)\n", "plt.axvline(CI_A[0], linestyle='--', color=\"blue\", ymin=0, ymax=1)\n", "plt.axvline(CI_A[1], linestyle='--', color=\"blue\", ymin=0, ymax=1)\n", "plt.axvline(CI_B[0], linestyle='--', color=\"orange\", ymin=0, ymax=1)\n", "plt.axvline(CI_B[1], linestyle='--', color=\"orange\", ymin=0, ymax=1)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another typical way to plot (only mean and CI):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot([1], mean_A, 'o', color=\"blue\")\n", "plt.plot([2], mean_B, 'o', color=\"orange\")\n", "plt.plot([1,1], CI_A, linestyle='-', color=\"blue\")\n", "plt.plot([2,2], CI_B, linestyle='-', color=\"orange\")\n", "plt.xlim([0,3])\n", "plt.xticks([1,2], [\"group A\", \"group B\"])\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Nutrition study 2:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff = mean_B - mean_A\n", "print(diff)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "se_diff = np.sqrt(se_A**2 + se_B**2)\n", "print(se_diff)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# Define the null hypothesis\n", "mean_null_hyp = 0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Compute the \"test statistic\" from the oberserved data\n", "tobs = (diff - mean_null_hyp) / se_diff\n", "print(tobs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# compute v (degrees of freedaom for difference-test)\n", "v = (se_A**2 + se_B**2)**2 / (se_A**4/(n_A-1) + se_B**4/(n_B-1))\n", "print(v)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# compare with t_0.975 from t-distribution with df = v\n", "stats.t.ppf(0.975, df=v)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "2*stats.t.cdf(-tobs, df=v)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# You can also use the ttest_ind function from scipy.stats:\n", "# \"_ind\" is for independent - i.e., two independent samples\n", "test = stats.ttest_ind(B,A,equal_var=False) \n", "# \"equal_var=False\" means we are NOT assuming equal variances - this is a \"Welch t-test\"\n", "print(test)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since 0.001 < p < 0.01 we conclude there is \"Strong evidence against H_0\"\n", "\n", "We reject the nullhypothesis\n", "\n", "And we conclude that the average energy usage of nurses in Hospital B is significantly larger than in Hospital A. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Nutrition study 3:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confidence interval:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diff_lower = diff - stats.t.ppf(0.975, df=v)*se_diff\n", "diff_upper = diff + stats.t.ppf(0.975, df=v)*se_diff\n", "print([diff_lower,diff_upper])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# You can also retrieve the confidence interval from the \"test\" calculated with stats.ttest_ind:\n", "test = stats.ttest_ind(B,A,equal_var=False) \n", "print(test.confidence_interval(0.95))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The confidence interval does not include the value 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Sleep medicine (paired t-test):" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "x1 = np.array([.7,-1.6,-.2,-1.2,-1,3.4,3.7,.8,0,2])\n", "x2 = np.array([1.9,.8,1.1,.1,-.1,4.4,5.5,1.6,4.6,3.4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(x1, alpha=0.5)\n", "plt.hist(x2, alpha=0.5)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "differences = x1-x2\n", "plt.hist(differences, alpha=0.5, color = \"green\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# make a \"one sample t-test\" for the differences (this time we only use the built-in python function):\n", "test1 = stats.ttest_1samp(differences, popmean=0)\n", "print(test1)\n", "print(test1.confidence_interval(0.95))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# equal to a paired t-test for two samples:\n", "test2 = stats.ttest_rel(x1, x2) # \"_rel\" is for related samples (= paired samples)\n", "print(test2)\n", "print(test2.confidence_interval(0.95))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# WRONG ANALYSIS:\n", "test3 = stats.ttest_ind(x1, x2) # making independent 2 sample test - WRONG!\n", "print(test3)\n", "print(test3.confidence_interval(0.95))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A paired experiment is \"stronger\" than an unpaired experiment. \n", "\n", "In a paired experiment each person is its owm \"negative control\" and we are sure to meassure a difference that is not due to inter-personal variation. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### QQ-plots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets re-visit the example with energy usage amongst nurses in hospital A and hospital B" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = np.array([7.53, 7.48, 8.08, 8.09, 10.15, 8.40, 10.88, 6.13, 7.90])\n", "B = np.array([9.21, 11.51, 12.79, 11.85, 9.97, 8.79, 9.69, 9.68, 9.19])\n", "\n", "plt.hist(A, alpha=0.5)\n", "plt.hist(B, alpha=0.5)\n", "plt.title(\"Hospital A (blue) and Hospital B (orange)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Could the data some from underlying normal distributions?\n", "\n", "Do we need to assume that the distributions are normal?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Lets make the QQ-plots:\n", "fig, axs = plt.subplots(1, 2, figsize=(10,4))\n", "sm.qqplot(A,line=\"q\",a=3/8,ax=axs[0]) \n", "sm.qqplot(B,line=\"q\",a=3/8,ax=axs[1]) \n", "# OBS: \"a = 3/8\" is preferred for n <= 10 \n", "# (\"a = 1/2\" is preferred for n > 10) \n", "axs[0].set_title(\"Hospital A\")\n", "axs[1].set_title(\"Hospital B\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Its difficult to tell whether it looks normal\n", "\n", "We can try to compare with simulated data from normal distributions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# for Hospital A\n", "fig, axs = plt.subplots(3, 3, figsize=(8,8))\n", "# data in first subplot:\n", "sm.qqplot(A,line=\"q\",a=3/8,ax=axs[0,0]) \n", "# simulated data in all other subplots:\n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[0,1]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[0,2]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[1,0]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[1,1]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[1,2]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[2,0]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[2,1]) \n", "sm.qqplot(stats.norm.rvs(size=len(A)),line=\"q\",a=3/8,ax=axs[2,2]) \n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# for Hospital B\n", "fig, axs = plt.subplots(3, 3, figsize=(8,8))\n", "# data in first subplot:\n", "sm.qqplot(B,line=\"q\",a=3/8,ax=axs[0,0]) \n", "# simulated data in all other subplots:\n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[0,1]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[0,2]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[1,0]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[1,1]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[1,2]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[2,0]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[2,1]) \n", "sm.qqplot(stats.norm.rvs(size=len(B)),line=\"q\",a=3/8,ax=axs[2,2]) \n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When comparing to simulated data, it looks like our data could very well come from an underlying normal distribution.\n", "\n", "This is no proof of normality(!), but it can give us an idea whether the assumption could be true." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Power calculations" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "# import package for power calculations:\n", "import statsmodels.stats.power as smp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# We want to calculate the sample size needed in a new experiment\n", "\n", "# from previous experiments we have a good guess for the sample variation:\n", "sd = 1.62\n", "\n", "# we want to be able to detect a voltage drop of down to 0.5 volts\n", "delta = 0.5\n", "\n", "# we need some fractals from the normal distribution (assuming the samples will be big enough for normal distribution assumption to apply)\n", "z_power = stats.norm.ppf(0.80, loc=0, scale = 1)\n", "z_signif = stats.norm.ppf(0.975, loc=0, scale = 1)\n", "\n", "# use formula to calculate n_obs (needed number of observations):\n", "n_obs = (sd/delta*(z_power+z_signif))**2\n", "print(n_obs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we can also use the python function TTestPower().solve_power (for one sample power calculations): \n", "print(smp.TTestPower().solve_power(effect_size=delta/sd, alpha=0.05, power=0.80))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results is almost the same (due to assumption anout normaldistribution)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# What if we can only make 60 observations?\n", "n = 60\n", "\n", "# calculate the power:\n", "z_power_new = np.sqrt(n*(delta/sd)**2)-z_signif\n", "print(z_power_new)\n", "\n", "power = stats.norm.cdf(z_power_new)\n", "print(power)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we can also use the python function TTestPower().solve_power (for one sample power calculations): \n", "print(smp.TTestPower().solve_power(effect_size=delta/sd, alpha=0.05, nobs=60))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is almost the same (due to assumption anout normaldistribution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example: Power calculations, 2 samples" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1)\n", "\n", "Find the sample size in a test where:\n", "\n", "power = 0.90 (beta = 0.10)\n", "\n", "n1 = n2 (k = 1)\n", "\n", "We want to be able to detect a difference of 2\n", "\n", "using alpha = 0.05\n", "\n", "and assuming sigma = 1 in both populations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "power = 0.90\n", "k = 1\n", "delta = 2\n", "sd = 1\n", "\n", "z_power = stats.norm.ppf(0.90, loc=0, scale = 1)\n", "z_signif = stats.norm.ppf(0.975, loc=0, scale = 1)\n", "\n", "n1 = (k+1) * (sd/delta*(z_signif+z_power))**2\n", "\n", "print(n1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we can also use the python function TTestIndPower().solve_power (for one sample power calculations): \n", "print(smp.TTestIndPower().solve_power(effect_size=delta/sd, alpha=0.05, power=0.90, ratio=k))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is almost the same (due to assumption anout normaldistribution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2)\n", "\n", "Find the power of an experiment where: \n", "\n", "n1 = n2 = 10 (k = 1)\n", "\n", "We want to be able to detect a difference of 2\n", "\n", "using alpha = 0.05\n", "\n", "and assuming sigma = 1 in both populations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n1 = 10\n", "n2 = 10\n", "k = 1\n", "delta = 2\n", "sd = 1\n", "z_signif = stats.norm.ppf(0.975, loc=0, scale = 1)\n", "\n", "z = np.sqrt(n1/(k+1)*delta**2/sd**2) - z_signif\n", "power = stats.norm.cdf(z)\n", "\n", "print(power)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we can also use the python function TTestIndPower().solve_power (for one sample power calculations): \n", "print(smp.TTestIndPower().solve_power(effect_size=delta/sd, alpha=0.05, nobs1=10, ratio=k))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is almost the same (due to assumption anout normaldistribution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3)\n", "\n", "What effectsize (delta) can be detected when:\n", "\n", "n1 = n2 = 10 (k = 1)\n", "\n", "alpha = 0.05\n", "\n", "power = 0.90 (beta = 0.10)\n", "\n", "and assuming sigma = 1 in both populations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n1 = 10\n", "n2 = 10\n", "k = 1\n", "sd = 1\n", "\n", "z_power = stats.norm.ppf(0.90, loc=0, scale = 1)\n", "z_signif = stats.norm.ppf(0.975, loc=0, scale = 1)\n", "\n", "delta = sd*(z_signif+z_power)/np.sqrt(n1/(k+1))\n", "\n", "print(delta)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we can also use the python function TTestIndPower().solve_power (for one sample power calculations): \n", "print(smp.TTestIndPower().solve_power(alpha=0.05, power=0.90, nobs1=10, ratio=k))" ] } ], "metadata": { "kernelspec": { "display_name": "pernille", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 2 }