--- title: 'Lecture 3: Estimates and confidence intervals' author: "Matti Pirinen" date: "11.8.2023" output: html_document: default --- {r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) set.seed(19)  Let's start the same way as in Lecture 2: Suppose we have 100 pairs of same-sex twins that are discordant for psoriasis and we know that BMI is higher for the psoriatic twin in 71 out of 100 pairs. Does that sound like higher BMI goes with psoriasis status in the population where these data come from? How to statistically quantify the association between BMI and psoriasis in these data? Previously, we quantified (or "tested") how consistent our observation was with a **null hypothesis** saying that "psoriasis status is independent of BMI". We did this by computing a P-value, that was the tail probability of observing at least as extreme an observation as what we had observed, assuming that the null hypothesis was true. This value gave an idea how clearly the observation is pointing to a possible deviation from the null hypothesis, with smaller P-values indicating more evidence against the null hypothesis. But P-value is only a statistical summary and typically a more interesting quantity is the actual value of the unknown variable that we want to estimate. Here that value is the proportion of twin pairs, where higher BMI and psoriasis go together. What can we say about the proportion of pairs where the psoriatic twin has the higher BMI? Should we say that it is 71%? Yes, in these data it happens to be exactly 71%, but what about more generally among a larger population of twin pairs. Maybe it would be 68% or 75%? One of the main goals of statistics is to make statements about the population (such as all twin pairs) based on a relatively small sample from that population (such as$n=100$twin pairs). Therefore, we need to quantify the uncertainty that exists because the sample we have observed is only a subset of the target population rather than the whole population. To demonstrate the sampling variation due to collecting only a subset of the population, let's simulate 10 data sets of 100 twin pairs under an assumption that the true proportion parameter in the target population is$p=0.71$. {r} p = 0.71 #true proportion in population n.trials = 100 #sample size in each data set rbinom(10, size = n.trials, prob = p) #returns 10 values, each btw 0,..,100.  We see that the observed value varies, and is exactly 71 in only one case above even if the true proportion in the simulation is exactly 0.71. To get a comprehensive picture of this statistical variation, let's generate 10,000 example data sets and show results as a table. (We don't want to print out vector of 10,000 values on screen but rather collect the result into a table where each observed value is represented only once.) {r} p = 0.71 #true proportion in population n.experiments = 10000 #no. of data sets n.trials = 100 #sample size in each dat set x = rbinom(n.experiments, size = n.trials, prob = p) #returns 10000 values, each btw 0,..,100. p.est = x / n.trials #10000 estimates ('est') of proportion parameter p table(p.est)/n.experiments #frequency distribution of estimates  In all these data sets, the true population parameter was 71%, but only in less than 1 in 10 of the data sets the observed proportion is exactly 71%. The observed proportion varies quite a bit due to random sampling (here range is from 0.52 to 0.88). Clearly, when we have observed a single **point estimate** of 71% from one data set of$n=100$samples, we need to report also something else than the point estimate 71%, to give an idea how much the true population proportion might vary from our point estimate due to sampling variation. Let's look at the distribution of point estimates when the true value is 0.71. {r} hist(p.est, breaks = 40, col = "limegreen") mean(p.est) sd(p.est)  Mean is where we would hope it to be, i.e., near the true value of 0.71. We say that our estimator is **unbiased**, that is, the estimates are "correct on average". Standard deviation of the estimates describes how much, on average, the estimate varies around its mean. We call this standard deviation as **standard error** (SE) of the estimator. That is, $$\text{SE of a parameter = "standard deviation of the sampling distribution of the parameter estimates".}$$ We also have a formula for the standard error for the mean of binomial sampling from Bin(n, p): $$\text{SE}_{\text{bin}} = \sqrt{\frac{p \,(1-p)}{n}}.$$ {r} SE = sqrt( p*(1-p) / n.trials) SE  Compare this to sd(p.est), which was an empirical estimate over 10000 samples. They do agree very well, so we are happy to use the formula in the future (except when$n$is small, when formula does not work that well). Let's see what is the interval within which, say, 95% of the point estimates fall, when the true parameter is 71%. Let's use quantile() function to get empirical cut-points of 2.5% and 97.5% of the distribution, since between them is 95% of the distribution. {r} quantile(p.est, c(0.025, 0.975) ) #we can give two cut-points the same time as a vector  So if the true parameter is 0.71, then in 19 out of 20 cases (that is 95% of cases) we see a point estimate in the range of 0.62,..,0.79. If we had not done this empirical simulation of 10000 binomial experiments, we could still have written down the standard 95% **confidence interval** (CI) estimate based on the SE, using the rule that 95% CI is given by (point.estimate +/- 1.96*SE). {r} c(0.71 - 1.96 * SE, 0.71 + 1.96 * SE) #Let's use 0.71 as our point estimate here.  This is a good approximation to 95% CI from the empirical data when n.trials is large (like in our case). Often it is enough to remember that 95% interval around the point estimate is given by adding 2 SEs on both sides of the estimate. More realistically, suppose that we do not know the true value for$p$, but we just have a point estimate p.est computed from the data. If we surround the estimate by its 95% CI computed from SE, what happens? {r} #currently p.est has 10000 point estimates, each from a data set of size 100 SE.est = sqrt( p.est * (1 - p.est) / n.trials) #vector of 10000 SE estimates, one for each data set ci.low = p.est - 1.96 * SE.est #vector of lower bounds of 95% CIs ci.up = p.est + 1.96 * SE.est #vector of upper bounds of 95% CIs  For example, for the first data set we have the following values computed. (Note that we can name values in a vector as follows.) {r} c(p.est = p.est[1], SE = SE.est[1], ci.low = ci.low[1], ci.up = ci.up[1])  Let's then visualize the first 100 point estimates with CIs, and compare them to the true value of 0.71 (red line). We first use plot(NULL) command to make an empty plot with suitable x-axis and y-axis limits (xlim and ylim parameters) and labels (xlab and ylab parameters). Then we plot all the point estimates using points() that takes in x-coordinates and y-coordinates and we use pch=19 to make solid plotting symbols and cex = 0.5 to make size 50% of normal. Finally, we draw CIs by using arrows() that takes in four vectors of coordinates (x-start, y-start, x-end, y-end) and draws lines between start and end points. code=0 means that lines have no symbols at ends and lwd = 0.5 means that line width is 50% of the normal. {r, fig.width = 10} n.to.plot = 100 #how many to plot -- plotting them all might be too messy plot(NULL, xlim = c(1, n.to.plot), ylim = c(0.4, 0.95), ylab = "p", xlab = "experiments") points(1:n.to.plot, p.est[1:n.to.plot], pch = 19, cex = 0.5) arrows(1:n.to.plot, ci.low[1:n.to.plot], 1:n.to.plot, ci.up[1:n.to.plot], code = 0, lwd = 0.5) abline(h = p, col = "red") #add horizontal line at y = p i.e. y=0.71  If we now would have observed only one of these 10000 data sets, and would report its point estimate and its 95% CI, in how many cases would the 95% CI cover the true value of 0.71? Let's make a logical TRUE / FALSE vector for each data set that is TRUE, if the 95% CI covers the true value, and is FALSE otherwise. {r} covers = (ci.low < p & ci.up > p) #TRUE if p is btw ci.low and ci.up, otherwise FALSE table(covers)/n.experiments  The coverage of 95% CIs here is quite close to 95% (actually slightly smaller here, about 93.6%). In general, when n.trials is large enough, about 95% of the CIs cover the true parameter value. This is the typical interpretation of the 95% confidence interval. In other words, in about 95% of the cases where we compute a point estimate and its 95% CI, the true parameter value is within the confidence interval. We cannot know for sure whether the true parameter value is within any one CI, but we know that, on average, it is there in 95% of the cases. We could also consider confidence intervals for other coverages than 95%. If we wanted a CI for coverage$\theta$, then we change the factor 1.96 multiplying SE to abs(qnorm((1-theta) / 2)). {r} theta = 0.95 #95% CI abs(qnorm((1-theta) / 2)) theta = 0.50 #50% CI abs(qnorm((1-theta) / 2))  Thus 50% confidence intervals are defined when we add 0.67*SE on both sides of the point estimates. Summary: A confidence interval around the point estimate gives an idea where the true population parameter may be. The interpretation of 95% CI is that if we were to repeat the experiment over and over again, then 95% of the CIs would cover the true parameter value and in the remaining 5% the CI would not cover the true value. When CIs are wide, we are very uncertain where the true parameter value may be and our experiment has simply not been very informative. It is crucial to report CIs in addition to the parameter estimate. Result "parameter is 5.6 (95%CI 0.0 .. 24.5)" is very different from "parameter is 5.6 (95%CI 5.2 .. 6.0)" even though both could be (uninformatively) summarized by "parameter estimate is 5.6". The former leaves us with a lot of uncertainty about the actual value of the parameter whereas the latter is more accurate. Often point estimate and SE or CI are more informative than the P-values under the null hypothesis (reading from [BMJ](https://www.bmj.com/content/346/bmj.f3212)), because estimates explicitly tell about the values of the parameters of interest whereas P-values tell only about statistical significance, which may not always have practical significance (as we saw [earlier](https://www.students4bestevidence.net/statistical-significance-vs-clinical-significance/)). #### Example 3.1 (1) Over the last year, there has been 316 appendectomy patients of which 180 were women and 136 were men. What is the point estimate of the proportion of women patients? What is the standard error of that estimate? How do you compute the 95% confidence interval for the estimate? How do you explain these results in plain language? Point estimate of the proportion of women is the observed number of women divided by the number of all patients. {r} n.all = 316 n.w = 180 p.w.est = n.w/n.all p.w.est #point estimate  SE comes from the formula above, where we replace the true proportion value with our current point estimate: {r} se = sqrt( p.w.est*(1-p.w.est)/n.all ) se  95% CI results when we add 1.96*SE on both sides of the point estimate: {r} c(p.w.est - 1.96*se, p.w.est + 1.96*se)  Explanation: We have observed that 180 out of 316 patients are women. We use these data to estimate the proportion of women among the whole patient population. The point estimate for the proportion of women patients is 180/316 = 57% (or equivalently 0.57). The 95% confidence interval is (51.5%,..,62.4%). This gives an idea of a range where the true population parameter might be based on the observed counts. The exact interpretation of the 95% confidence interval is that in 95% of the data sets, the CI covers the true parameter value but in the remaining 5% the true value lies outside the CI. (2) How do the values from (1) change if all sample sizes are multiplied by ten, i.e., there are 3160 patients and 1800 women and 1360 men? Explain intuitively why something changes and why something else does not change. {r} n.all = 3160 n.w = 1800 p.w.est = n.w/n.all p.w.est #point estimate se = sqrt( p.w.est*(1-p.w.est)/n.all ) se c(p.w.est - 1.96*se, p.w.est + 1.96*se)  The point estimates stays exactly the same when counts are multiplied by 10. Standard error gets smaller, reflecting the fact that large sample size$n$means less error in estimate. Similarly, 95%CI gets narrower reflecting the fact that with larger sample we have more information and therefore more accurate estimates of the parameters. #### Confidence intervals in binom.test() Previously, we did the binomial test using the function binom.test(). Let's revisit it now to see what it says about confidence intervals. {r} binom.test(71, n = 100, p = 0.5)  In addition to the P-value, it also gives 95% CI and the point estimate ("sample estimate"). We could alter the coverage of the CI by parameter conf.level that, by default, is 0.95. {r} binom.test(71, n = 100, p = 0.5, conf.level = 0.5)  Next we move to more general tests about proportion parameters. ### Proportion test We have used binom.test() to do binomial test. Another and more general function is prop.test() that is only an approximation to the binomial distribution but can be used also to compare two or more samples against each other. Let's study prop.test() in more detail. Type ?prop.test on console to read help. **prop.test for one group.** Let's first revisit question of 71 successes out of 100 trials. First with the two-sided test and then with the one sided test that assumes that the proportion is > 0.5. {r} prop.test(71, n = 100, p = 0.5, alternative = "two.sided", conf.level = 0.95) #And with one-sided test prop.test(71, n = 100, p = 0.5, alternative = "greater", conf.level = 0.95)  In the output, X-squared is the value of test statistics that leads to the P-value and df is "degrees of freedom" of the test (here 1 as there is one free parameter, p). The P-values are close to the exact values given by binom.test(). Note how null value is the same, 0.5, in both tests but the alternative hypothesis is different. As expected, 2-sided P-value is twice the 1-sided P-value. Note also how the confidence interval changes with the alternative hypothesis. (It is not recommended to use CIs for one-sided alternatives as their interpretation is complex. As a general rule, use two-sided tests and CIs.) **prop.test for two groups.** The benefit of prop.test() over binom.test() is that we can also compare two samples against each other. For example, suppose that in a similar study design in Norway it has been observed that the twin with psoriasis has higher BMI in 344 times out of 497 pairs of twins. Let's make a data matrix, where rows are countries and column 1 counts when psoriasis goes with higher BMI and column 2 the opposite cases. {r} x = rbind(FIN = c(71, 100-71), NOR = c(344, 497-344)) #bind rows together x #show the data to make sure we got it correct prop.test(x)  prop.test() runs with such a 2x2 matrix and tests whether the proportions in rows are different from each other. In this case, there is no statistically significant difference between the propotions in Finland (0.710) and Norway (0.692). The CI given in output is for the difference in the two proportions. So we may say that "the point estimate for the difference between Finland and Norway is 0.71-0.692 = 0.018 and its 95% confidence interval is from -0.086 to 0.122. In particular, 0 is pretty much in the middle of the CI and hence there is no evidence that the difference would be different from 0. The P-value also tells that there is no evidence for deviation from the null hypothesis that the proportions are the same in the two populations (P = 0.8144). ----- **Difference between proportions.** (This is to explain what prop.test() does under the hood, but this is not needed in the exercises.) Consider two populations and assume that proportions of cases (e.g. carriers of a particular disease) are$p_1$and$p_2$. We have access to the samples of sizes$n_1$and$n_2$and estimate$\widehat{p}_1 = x_1/n_1$and$\widehat{p}_2 = x_2/n_2$where$x_i$is the number of observed cases in population$i=1,2.$("Hat"-notation on top of$p_1$and$p_2$means the estimates that we have made about these parameters. We do not know the true values of$p_1$and$p_2$but we can estimate them using observed data by$\widehat{p}_1$and$\widehat{p}_2$defined above.) We compute the difference$\widehat{d} = \widehat{p}_1 - \widehat{p}_2$and want to do statistical inference about that difference. We have that $$\textrm{SE}\left(\widehat{p}_1 - \widehat{p}_2\right) \approx \sqrt{ \frac{\widehat{p}_1(1-\widehat{p}_1)}{n_1} + \frac{\widehat{p}_2(1-\widehat{p}_2)}{n_2}},$$ For testing the null hypothesis that$p_1=p_2$, we use the test statistic $$\frac{\widehat{p}_1 - \widehat{p}_2 }{\sqrt{\widehat{p}(1-\widehat{p})\left(\frac{1}{n_1}+\frac{1}{n_2} \right)}}, \textrm{ where } \widehat{p} = \frac{x_1+x_2}{n_1+n_2}.$$ Under the null hypothesis$p_1=p_2$, this test statistic follows approximately a standard Normal distribution (Topic of the next lecture), from which we can derive a P-value. ----- #### Example 3.2 (1) Over the last year, in Hospital A, there has been 316 appendectomy patients of which 180 women and 136 men. Use prop.test() to test whether both sexes are equally represented among appendectomy patients. {r} #Here we run a one-sample prop.test by specifying the observed number of women, # total sample size 'n' and null hypothesis success rate 'p' prop.test(180, n = 316, p = 0.5)  (2) In another hospital B, there has been 245 patients of which 134 women and 111 men. Use prop.test() to test whether the proportion of women is different between the two hospitals A and B. {r} #We make a matrix where one row is one hospital and columns are # numbers of women and men. x = rbind(A = c(180, 136), B = c(134, 111)) colnames(x) = c("women","men") x prop.test(x)  We see that the observed proportions of women are 0.570 and 0.547, which are not statistically different from each other (P-value 0.65) and the 95%CI for the difference is (-0.064, 0.109), containing 0 within it. **prop.test for more than two groups.** What if we have more than two groups to compare? Suppose three Finnish hospitals have done knee replacement operations in 2013 and in two year check-up the patients reported knee pain as follows. | Hospital | no pain | pain | total | |:--------:|:-------:|:----:|:-----:| |A|113|198|312| |B|100|111|211| |C|207|215|422| Are there differences between the hospitals? {r} x = rbind(A = c(113,198), B = c(100,111), C = c(207,215)) colnames(x) = c("nopain","pain") x prop.test(x)  Here prop.test() did a comparison to null hypothesis that all proportions are the same. A small P-value alone is not very informative since it does not tell which of the proportions look different from other(s). Therefore, pairwise comparisons might be more informative here. Note that here we can't estimate a single value for the "difference between groups" since there are more than 2 groups and hence we don't have a confidence interval either. ### Chi-square test Suppose we divide painless patients into two groups: fully functional and partially functional making a 3x3 table | Hospital | func+ | func- | pain | total | |:--------:|:-----:|:-----:|:----:|:-----:| |A|62|51|198|312| |B|50|50|111|211| |C|84|123|215|422| {r} x = rbind(A = c(62,51,198), B = c(50,50,111), C = c(84,123,215)) colnames(x) = c("func+", "func-", "pain") x  We can't use prop.test() anymore since there are more than two columns, that is, we are not observing a single proportion per hospital anymore. Now we should ask more generally whether every row has the same distribution, that is, are the proportions of the outcomes same across hospitals. Another way to put this is to ask whether the rows and columns are *independet*, that is, whether the distribution of counts in a cell on row$i$column$j$in the table is determined by the product of probabilities of row$i$and of column$j$, for all$i$and$j$. We can test this by a general **chi-square test** for contingency tables using chisq.test(). Read ?chisq.test. {r} chisq.test(x)  chisq.test() takes in a matrix and returns a P-value under the null hypothesis that the rows and columns are independent of each other. Here it is saying that the rows and columns do not seem independent (P-value is small), but otherwise it is not very informative, since it does not indicate where the possible differences are. More useful may be at least to see the numerical proportions per hospital to consider which groups and in which hospitals seem to differ from the other hospitals. {r} x/rowSums(x)  From these proportions it looks like the hospital A may be different from B and C. ----- **How chisq.test() does the chi-square ($\chi^2$) test for a contingency table?** (1) Calculate the expected frequency ($E_{ij}$) for the observation in row$i$and column$j$of the$r \times c$contingency table:$E_{ij} = \frac{i\textrm{th row total}\, \times\, j\textrm{th col total}}{\textrm{table total}}$(2) For each cell in the table calculate the difference between the observed value and the expected value$(O_{ij} - E_{ij})$. (3) Square each difference and divide the resultant quantity by the expected value$(O_{ij}-E_{ij})^2/E_{ij}$. (4) Sum all of these$r\times c$values to get a single number, the$\chi^2$statistic X.sq. (5) Compare this number with the chi-squared distribution with the following degrees of freedom:$\textrm{df} = (r - 1) \times (c - 1)\$. P-value is in R given by 1-pchisq(X.sq, df = df). Chi-square test is applicable if almost all expected cell counts are at least 5. If there are smaller counts, then exact methods such as Fisher's test should be used instead. ----- #### Example 3.3 (1) Over the last year, in Hospital A, there has been 316 appendectomy patients of which 180 women and 136 men. In another hospital B, there has been 245 patients of which 134 women and 111 men. Use Chi-square test to test independence between rows (hospitals) and columns (sex). Compare the P-value to one from prop.test(). {r} x = rbind(A = c(180,136), B = c(134,111)) colnames(x) = c("women","men") x chisq.test(x) prop.test(x)  Thus, in case where both prop.test() and chisq.test() are possible, they give the same test result. But prop.test() is more informative since it also tells the CI for the difference between the two proportions as well as the point estimates of the two proportions. **Summary of the tests we have used so far:** (1) binom.test() does a single sample binomial test, that is, compares the observed number of successes to the expectation under the null hypothesis success probability value, and returns a P-value, a point estimate and a confidence interval for the proportion. (2) prop.test() is an approximation to binom.test() in case of one sample, but can also compare many proportions against each other, not only to the null value as binom.test(). For the case of two groups, it gives a CI for the difference between the two success proportions. (3) chisq.test() takes in any matrix and compares it to the null hypothesis that all the rows and columns are independent of each other, that is, all the rows have the same distribution of how the observations are distributed into the columns. A very general test, but returns only a P-value and not any informative parameter about what kind of a difference from the null hypothesis might have been observed. (4) fisher.test() is suitable for tables that have small counts as there it is exact rather than an approximation as chisq.test(). In tables with larger counts, it gives similar results to chisq.test().