Dear R-ussers,

While looking for a way to calculate the association between a countinuous 
and a binary variable, 
i found a procedure called point biserial corralation.
Me, not being a mathematicion, i did my very best to understand what it 
was all about, and then i found a easily understandable paper (by steve 
simon) on ow to calculate this. ref ## 
http://www.childrens-mercy.org/stats/definitions/biserial.htm (this page 
has the same example)
Further i discovered the polycor package in R.
Now i'm having troubles with the fact that the polycor pkg never gives me 
the same output as the manuals aplication of the formula.

In the example below found, manualy  r(biserial) = 0.49 between fb an age, 
and ussing function polyserial (polycor pkg)  r(biserial) =-0.8591.
This is a rather big difference, no due to abriviation or flootingpoints.

Is there someone whom is familiar with biserial correlation, and the 
appropriate way to calculate it?

Kind regards,
Tom.

here is the example, at the end is the R file. 

1e I create the input

> library(abind) 
> library (polycor)
> ### data input
> no <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8)
> fb <- c(19, 30, 20, 19, 29, 25, 21, 24, 50, 25, 21, 17, 15, 14, 14, 22, 
17) 
> ss <- c(14, 41, 18, 11, 16, 24, 18, 21, 37, 17, 10, 16, 22, 12, 14, 12, 
18)
> age <-c('elderly', 'elderly', 'elderly', 'elderly', 'elderly', 
'elderly', 'elderly', 'elderly', 'elderly', 'young', 'young', 'young', 
'young', 'young', 'young', 'young', 'young') 
> 
> dataset <- data.frame(no,fb,ss,age)
> dataset <- subset(dataset,select=c(fb:age))
> nrow(dataset) 
[1] 17
> data_eld <- subset(dataset,age=='elderly',select=fb)
> data_young <- subset(dataset,age=='young',select=fb)
> 

here i calculate the R_bis (biserial corelation) manualy

> ### point biserial correlation
> 
>  fb <- subset(dataset,select=fb)
>  fb0 <- subset(dataset,age!='elderly',select=fb)
>  fb1 <- subset(dataset,age=='elderly',select=fb)
>  meanfb0 <- mean(fb0,na.rm=T)
>  meanfb1 <- mean(fb1,na.rm=T)
>  sdfb<- sd(dataset$fb,na.rm=T)
> 
>  ss <- subset(dataset, select=ss)
>  ss0 <- subset(dataset,age!='elderly',select=ss)
>  ss1 <- subset(dataset,age=='elderly',select=ss)
>  meanss0 <- mean(ss0,na.rm=T)
>  meanss1 <- mean(ss1,na.rm=T)
>  sdss<- sd(dataset$ss,na.rm=T) 
>  age <- subset(dataset,select=age) 
>   n <- nrow(dataset)
> 
this is the formula from  ref ## 
http://www.childrens-mercy.org/stats/definitions/biserial.htm

> > R_bis <- function(x,x1,x0,n){  p <- (nrow(x1)/n)
>+        (((mean(x1,na.rm=T)-mean(x0,na.rm=T))/sd(x,na.rm=T)) 
*sqrt(p*(1-p)))  } 

this is the corrected formula from  ref ## 
http://en.wikipedia.org/wiki/Point-biserial_correlation_coefficient
>> R_bis2 <- function(x,x1,x0,n){ 
+        ((mean(x1,na.rm=T)-mean(x0,na.rm=T))/sd(x,na.rm=T))  * (sqrt( 
(nrow(x1)*nrow(x0))/(n*(n-1))))}
> 
>>  R_bis(fb,fb1,fb0,n) 
>       fb 
>0.4798873 

result in paper was 0.49 
 
>>   R_bis2(fb,fb1,fb0,n)
>       fb 
>0.4946565 

equals result in paper  0.49 

Then i use the polycor package,
function hetcor will give all the different correlation ressults

 
>> hetcor(dataset$fb,dataset$ss,dataset$age ,ML=TRUE)
>
>Maximum-Likelihood Estimates
>
>Correlations/Type of Correlation:
 >           dataset$fb dataset.ss dataset.age
>dataset$fb           1    Pearson  Polyserial
>dataset.ss       0.703          1  Polyserial
>dataset.age    -0.8591    -0.6685           1
>
>Standard Errors:
>           dataset$fb dataset.ss
>dataset$fb 
>dataset.ss      0.1215 
>dataset.age     0.1106     0.2497
>
>n = 17 
>
>P-values for Tests of Bivariate Normality:
>            dataset$fb dataset.ss
>dataset$fb 
>dataset.ss      0.1782 
>dataset.age     0.4269     0.4034
> hetcor(dataset,ML=TRUE)
>
>Maximum-Likelihood Estimates
>
>Correlations/Type of Correlation:
 >        fb      ss        age
>fb        1 Pearson Polyserial
>ss    0.703       1 Polyserial
>age -0.8591 -0.6685          1
>
>Standard Errors:
 >       fb     ss
>fb 
>ss  0.1215 
>age 0.1106 0.2497
>
>n = 17 
>
>P-values for Tests of Bivariate Normality:
>       fb     ss
>fb 
>ss  0.1782 
>age 0.4269 0.4034

here a quick two step method is ussed to calculate the polyserial 
correlation

> >polyserial(dataset$fb,dataset$age)
>[1] -0.6205737
>> polyserial(dataset$fb,dataset$age, ML=TRUE, std.err=TRUE) 

same method  as in hetcor, only for indecated variables

>Polyserial Correlation, ML est. = -0.8591 (0.1106)
>Test of bivariate normality: Chisquare = 4.91, df = 5, p = 0.4269
>
>               1
>Threshold 0.1811
>Std.Err.  0.1849
> 


> ### for side to side (ss)   incase no 9 is an outlier in fb, this will 
not be the case in ss 
> 
>> R_bis(ss,ss1,ss0,n)
>       ss 
>0.4153681 

result in paper was 0.43 

>> R_bis2(ss,ss1,ss0,n)
>       ss 
>0.4281516 

equals result in paper  0.43

>> polyserial(dataset$ss,dataset$age) 
>[1] -0.5371397

>> polyserial(dataset$ss,dataset$age, ML=TRUE, std.err=TRUE)
>
>Polyserial Correlation, ML est. = -0.6685 (0.2497)
>Test of bivariate normality: Chisquare = 5.103, df = 5, p = 0.4034
>
>               1
>Threshold 0.1504
>Std.Err.  0.2583
#################################################################################

### R-file
#################################################################################

library(abind) 
library (polycor)
### data input
no <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8)
fb <- c(19, 30, 20, 19, 29, 25, 21, 24, 50, 25, 21, 17, 15, 14, 14, 22, 
17) 
ss <- c(14, 41, 18, 11, 16, 24, 18, 21, 37, 17, 10, 16, 22, 12, 14, 12, 
18)
age <-c('elderly', 'elderly', 'elderly', 'elderly', 'elderly', 'elderly', 
'elderly', 'elderly', 'elderly', 'young', 'young', 'young', 'young', 
'young', 'young', 'young', 'young') 

dataset <- data.frame(no,fb,ss,age)
dataset <- subset(dataset,select=c(fb:age))
nrow(dataset) 
data_eld <- subset(dataset,age=='elderly',select=fb)
data_young <- subset(dataset,age=='young',select=fb)

### point biserial correlation

 fb <- subset(dataset,select=fb)
 fb0 <- subset(dataset,age!='elderly',select=fb)
 fb1 <- subset(dataset,age=='elderly',select=fb)
 meanfb0 <- mean(fb0,na.rm=T)
 meanfb1 <- mean(fb1,na.rm=T)
 sdfb<- sd(dataset$fb,na.rm=T)

 ss <- subset(dataset, select=ss)
 ss0 <- subset(dataset,age!='elderly',select=ss)
 ss1 <- subset(dataset,age=='elderly',select=ss)
 meanss0 <- mean(ss0,na.rm=T)
 meanss1 <- mean(ss1,na.rm=T)
 sdss<- sd(dataset$ss,na.rm=T) 
 age <- subset(dataset,select=age) 
  n <- nrow(dataset)
 
 R_bis <- function(x,x1,x0,n){  p <- (nrow(x1)/n)
       (((mean(x1,na.rm=T)-mean(x0,na.rm=T))/sd(x,na.rm=T)) 
*sqrt(p*(1-p)))  }
R_bis2 <- function(x,x1,x0,n){ 
       ((mean(x1,na.rm=T)-mean(x0,na.rm=T))/sd(x,na.rm=T))  * (sqrt( 
(nrow(x1)*nrow(x0))/(n*(n-1))))}
 
 R_bis(fb,fb1,fb0,n)
  R_bis2(fb,fb1,fb0,n)

 
hetcor(dataset$fb,dataset$ss,dataset$age ,ML=TRUE)
hetcor(dataset,ML=TRUE)

polyserial(dataset$fb,dataset$age)
polyserial(dataset$fb,dataset$age, ML=TRUE, std.err=TRUE) 

### for side to side (ss)   incase no 9 is an outlier in fb, this will not 
be the case in ss 

R_bis(ss,ss1,ss0,n)
R_bis2(ss,ss1,ss0,n)

polyserial(dataset$ss,dataset$age) 
polyserial(dataset$ss,dataset$age, ML=TRUE, std.err=TRUE)

## END 
##########################################################################

Willems Tom

E-mail: towil@var.fgov.be




Disclaimer: click here
	[[alternative HTML version deleted]]

