Chapter 2-14. Pearson Correlation Coefficient With Clustered Data
The ordinary Pearson correlation coefficient, as well as the nonparametric correlation coefficients, assume that each x-y pair of numbers is independent from the other x-y pairs. Both the interpretation of the correlation coefficient itself, as well as the significance test for the coefficient, require this assumption.
Bland and Altman published three one-page articles discussing this situation.
In their first paper (Bland and Altman, 1994), they simply present the idea that it is incorrect to compute a correlation coefficient on data that uses more than one x-y pair per study subject (the repeated measurement situation).
In their second paper (Bland and Altman, 1995a), they describe how to analyze such data if the interest is a “within subjects” research question. For example, suppose you measured pH and Paco2 in the same patient, and you repeated this approximately four times during their hospital stay. You have a “within subjects” research question if you want to know whether an increase in pH within the same individual was associated with an increase in Paco2. (Paco2 is the symbol for partial pressure of carbon dioxide in the arterial blood.)
In their third paper (Bland and Altman, 1995b), they describe how to analyze such data if the interest is a “between subjects” research question. For this same dataset, you have a “between subjects” research question if you want to know whether subjects with high values of pH also tend to have high values of Paco2.
In this chapter, both Bland and Altman’s within subjects approach and their between subjects approach is described, along with some Stata programs for doing them.
______
Source: Stoddard GJ. Biostatistics and Epidemiology Using Stata: A Course Manual [unpublished manuscript] University of Utah
School of Medicine, 2011. http://www.ccts.utah.edu/biostats/?pageId=5385
Significance test for Pearson Correlation Coefficient
We will need some Stata code for computing the significance test for the correlation coefficient when we get to the “between subjects” analysis below. Stata does not have a command where you can provide Stata with a correlation coefficent and sample size, after which it gives you the p value.
The hypothesis test that the population correlation coefficient is different from zero is (Rosner, 2006, p.496):
H0: ρ = 0 vs H1: ρ ¹ 0 , where ρ is the population correlation coefficient
test statistic: , where r is the sample correlation coefficient
For a two-sided level α test, reject H0 if
, where n is the sample size
Thus, for a two-sided α = 0.05 comparison, statistical significance (p < 0.05) is achieved when
To compute the two-tailed, or two-sided comparison, p value for a Pearson correlation coefficient, for a given sample size, cut-and-paste the following into the Stata do-file editor and run it. This loads the program to run in your current Stata session.
* program to compute 2-tailed p value, passing* it the Pearson r and the sample size
capture program drop pearsonpval
program define pearsonpval
args r n
local t=`r'*sqrt(`n'-2)/sqrt(1-`r'^2)
local p=2*ttail(`n'-2,abs(`t'))
disp as result _newline "r = " `r'
disp as result "n = " `n'
disp as result "p = " %6.4f `p'
end
* syntax for use: pearsonpval r n
Now, to compute the p value for an r = 0.25 based on a sample size of n = 100, you use
r = .25
n = 100
p = 0.0121
This result agrees with Rosner (2006, p.497), so you can feel confident it is programmed correctly.
Within Subjects Pearson Correlation Coefficient for Repeated Measures Data
As stated above, Bland and Altman (1995a) describe how to analyze repeated measures data if the interest is a “within subjects” research question. For example, you have a “within subjects” research question if you want to know whether an increase in pH within the same individual was associated with an increase in Paco2.
This section will work through the process described in Bland and Altman (1995a).
Bland and Altman (1995a) illustrate with a dataset they provide in their article (their Table 1). To read this dataset into Stata,
FileOpen
Find the directory where you copied the course CD:
Change to the subdirectory datasets & do-files
Single click on bland&altmanBMJ1995Table1.dta
Open
use "C:\Documents and Settings\u0032770.SRVR\Desktop\
regressionclass\datasets & do-files\
bland&altmanBMJ1995Table1.dta", clear
* which must be all on one line, or use:
cd "C:\Documents and ettings\u0032770.SRVR\Desktop"
cd "\regressionclass\datasets & do-files"
use bland&altmanBMJ1995Table1, clear
Listing the first three of the n=8 subjects,
list , sepby(subject) , if subject<=3+------+
| subject ph paco2 |
|------|
1. | 1 6.68 3.97 |
2. | 1 6.53 4.12 |
3. | 1 6.43 4.09 |
4. | 1 6.33 3.97 |
|------|
5. | 2 6.85 5.27 |
6. | 2 7.06 5.37 |
7. | 2 7.13 5.41 |
8. | 2 7.17 5.44 |
|------|
9. | 3 7.4 5.67 |
10. | 3 7.42 3.64 |
11. | 3 7.41 4.32 |
12. | 3 7.37 4.73 |
13. | 3 7.34 4.96 |
14. | 3 7.35 5.04 |
15. | 3 7.28 5.22 |
16. | 3 7.3 4.82 |
17. | 3 7.34 5.07 |
+------+
We see that each subject has multiple observations of pH and Paco2, with a variable number of observations per subject.
Suppose we are interested in a “within subjects” research question, where we want to know whether an increase in pH within the same individual was associated with an increase in Paco2.
This can be done with analysis of covariance (ANCOVA). We make one of the two variables the outcome variable and the other the continuous predictor variable. In addition, we make subject a predictor variable, by using indicator variables (7 indictor, or dummy, variables to represent the 8 subjects). The anova command automatically creates these indictor variables, without showing them to you, for any variable that is not listed in the continuous( ) option.
anova ph paco2 subject , continuous(paco2)Number of obs = 47 R-squared = 0.8993
Root MSE = .093715 Adj R-squared = 0.8781
Source | Partial SS df MS F Prob > F
------+------
Model | 2.98016151 8 .372520188 42.42 0.0000
|
paco2 | .115324822 1 .115324822 13.13 0.0008
subject | 2.96606639 7 .42372377 48.25 0.0000
|
Residual | .333732534 38 .008782435
------+------
Total | 3.31389404 46 .072041175
This ANCOVA table shows how the variability in pH can be partitioned into components due to different sources.
To show graphically what the model did, we ask for predicted values from the model
and then overlay line graphs through these predicted values for each subject onto the scatterplot of observations,
#delimit ;twoway (scatter ph paco2 , symbol(circle) color(black))
(line ph_pred paco2 if subject==1, color(black) clstyle(solid))
(line ph_pred paco2 if subject==2, color(black) clstyle(solid))
(line ph_pred paco2 if subject==3, color(black) clstyle(solid))
(line ph_pred paco2 if subject==4, color(black) clstyle(solid))
(line ph_pred paco2 if subject==5, color(black) clstyle(solid))
(line ph_pred paco2 if subject==6, color(black) clstyle(solid))
(line ph_pred paco2 if subject==7, color(black) clstyle(solid))
(line ph_pred paco2 if subject==8, color(black) clstyle(solid))
, legend(off) ytitle("Intramural pH") xtitle(PaCO2)
;
#delimit cr
We see that the model fit parallel lines through each subject’s data.
The residual sum of squares in the ANCOVA table (Residual SS, 0.333732534) represents the variation around these regression lines. The between subjects variation, along with any unmeasured “nuisance variables” which make subjects differ from each other, is reflected in the subjects term of the ANCOVA table (Subjects SS, 2.96606639). We remove that source of variation by not including it in our estimation of the within-subjects correlation.
We express the within subjects variation in observations of pH due to Paco2 as a proportion of what’s left:
This term is analogous to a coefficient of determination. When you square a Pearson correlation coefficient, you have the “coefficient of determination (r2)” which represents the proportion of variation in the outcome variable that is explained by the predictor variable.
To obtained the magnitude of the “correlation coefficient within subjects”, we take the square root of this proportion.
Filling in the equation from the ANCOVA table,
Number of obs = 47 R-squared = 0.8993
Root MSE = .093715 Adj R-squared = 0.8781
Source | Partial SS df MS F Prob > F
------+------
Model | 2.98016151 8 .372520188 42.42 0.0000
|
paco2 | .115324822 1 .115324822 13.13 0.0008
subject | 2.96606639 7 .42372377 48.25 0.0000
|
Residual | .333732534 38 .008782435
------+------
Total | 3.31389404 46 .072041175
The sign of the correlation coefficient comes from the sign of the regression coefficient for Paco2 when this ANCOVA is obtained from linear regression. This model fitted with linear regression is
i.subject _Isubject_1-8 (naturally coded; _Isubject_1 omitted)
Source | SS df MS Number of obs = 47
------+------F( 8, 38) = 42.42
Model | 2.98016151 8 .372520188 Prob > F = 0.0000
Residual | .333732534 38 .008782435 R-squared = 0.8993
------+------Adj R-squared = 0.8781
Total | 3.31389404 46 .072041175 Root MSE = .09371
------
ph | Coef. Std. Err. t P>|t| [95% Conf. Interval]
------+------
paco2 | -.108323 .0298928 -3.62 0.001 -.1688378 -.0478082
_Isubject_2 | .7046112 .0773549 9.11 0.000 .5480145 .861208
_Isubject_3 | .9500127 .0610954 15.55 0.000 .8263315 1.073694
_Isubject_4 | .9715577 .0735091 13.22 0.000 .8227464 1.120369
_Isubject_5 | .8603818 .0583954 14.73 0.000 .7421664 .9785972
_Isubject_6 | .9264284 .0659945 14.04 0.000 .7928296 1.060027
_Isubject_7 | .6921055 .1049093 6.60 0.000 .4797277 .9044834
_Isubject_8 | .7033361 .0615714 11.42 0.000 .5786913 .8279809
_cons | 6.929854 .129469 53.53 0.000 6.667758 7.19195
------
From the coefficient for paco2 (B = -.108323), we see the assocation is negative.
The p value for the correlation coefficient can be taken from the t test for the regression slope (p = 0.001) or from the F test for the paco2 term in the ANCOVA table ( p = 0.0008).
Number of obs = 47 R-squared = 0.8993
Root MSE = .093715 Adj R-squared = 0.8781
Source | Partial SS df MS F Prob > F
------+------
Model | 2.98016151 8 .372520188 42.42 0.0000
|
paco2 | .115324822 1 .115324822 13.13 0.0008
subject | 2.96606639 7 .42372377 48.25 0.0000
|
Residual | .333732534 38 .008782435
------+------
Total | 3.31389404 46 .072041175
It makes no difference which variable we choose for the dependent variable, as either way gives an identical correlation coefficient and p value.
These steps are contained in the following program. Simply cut-and-paste this into the do-file editor, and run it. This sets up the “withincorr” commmand to work during your current Stata session.
* Bland and Altman (1995)
* correlation coefficient within subjects
capture program drop withincorr
program define withincorr
version 10
syntax varlist(min=3 max=3) [if] [in]
tokenize `varlist'
local yvar `1'
local xvar `2'
local subjectid `3'
tempname touse
mark `touse' `wgt' `if' `in'
preserve
quietly keep if `touse'
quietly anova `yvar' `xvar' `subjectid' ///
, continuous(`xvar')
local r = sqrt(e(ss_1)/(e(ss_1)+e(rss)))
quietly xi: regress `yvar' `xvar' i.`subjectid'
if _b[`xvar']<0 {
local r = -`r'
}
local p = 2*ttail(`e(N)'-2,abs(_b[`xvar']/_se[`xvar']))
display as result _newline "within r = " ///
%4.2f `r' " , n = " `e(N)' " , p = " %5.4f `p'
restore
end
* syntax: withincorr y x subjectid
* ------
After executing the above lines to set up the “withincorr” command for the duration of your current Stata session, to compute a within subjects correlation coefficient, you run the following command,
Syntax: withincorr yvar xvar subjectid [the subject ID variable must come last]
withincorr ph paco2 subjectwithin r = -0.51 , n = 47 , p = 0.0007
This is what is found in the Bland and Altman (1995a) paper, so you can feel confident it was programmed correctly.
The command also works with the “if” option and “in” option. To illustrate,
withincorr ph paco2 subject if subject<5withincorr ph paco2 subject in 1/10
within r = -0.07 , n = 22 , p = 0.7628
within r = 0.02 , n = 10 , p = 0.9656
Here is a suggested way to report this in your article.
Article Suggestion
Statistical Methods
Pearson correlation coefficients are reported. For variables involving repeated measurements for each study subject, an ordinary correlation coefficient is not appropriate (Bland and Altman, 1994). For these variables, a “within subjects” correlation coefficient was reported, which accounts for the lack of independence among the repeated measurements by removing the variation between subjects (Bland and Altman, 1995a). A within subjects correlation coefficient examines whether an increase in a variable within the same individual is associated with an increase in another variable (Bland and Altman, 1995a).
Results
In the Results section, just report it like you would any other correlation coefficient. I doubt the reader would want to be bothered with knowing which correlation coefficients are ordinary coefficients and which are within subjects coefficients when reading the results.
Between Subjects Pearson Correlation Coefficient for Repeated Measures Data
In their third paper (Bland and Altman, 1995b), the authors describe how to obtain a correlation coefficient for repeated measures data if the interest is a “between subjects” research question. Using the same dataset provided for ther within subjects correlation, you have a “between subjects” research question if you want to know whether subjects with high values of pH also tend to have high values of Paco2.
Bland and Altman (1995b) illustrate with a dataset they provided in their second article (Bland and Altman, 1995a).
To read this dataset into Stata,
FileOpen
Find the directory where you copied the course CD:
Change to the subdirectory datasets & do-files
Single click on bland&altmanBMJ1995Table1.dta
Open
use "C:\Documents and Settings\u0032770.SRVR\Desktop\
regressionclass\datasets & do-files\
bland&altmanBMJ1995Table1.dta", clear
* which must be all on one line, or use:
cd "C:\Documents and ettings\u0032770.SRVR\Desktop"
cd "\regressionclass\datasets & do-files"
use bland&altmanBMJ1995Table1, clear
Listing the first three of the n=8 subjects,
list , sepby(subject) , if subject<=3+------+