y <- c(26073,18103,15100,14738,22259,10277,6225,7663,6618,9945,7590,4293,4656,5145,2001,1779,1609,2191,1617,2554,2198,1578,3446,1380,1249,1223,834,3754,2283,3028,1100,457,1201,2192,1508,1393,952,1032,1279,1370,649,1900,666,1313,1353,1500,877,874,1133,754,695,609,696,756,670,301,630,798,436,388,864,497,449,919,536,673,837,534,845,626,871,740,391,435,424,338,744,368,393,938,804,456,267,338,NA,635,402,462,525,1069,487,756,360,481,NA,480,68,262,520,298) x <- c(22274,14819,15136,13704,19638,7551,8019,6509,6634,11166,7508,4275,4944,5441,1689,1522,1416,1594,1909,2599,1262,1199,4404,1166,1122,886,778,4436,1890,3107,1038,300,988,2008,1522,1336,976,798,869,1260,578,2359,736,1690,1201,813,778,687,1270,671,1559,489,773,629,637,277,776,1651,377,222,1068,399,547,668,451,724,853,434,730,612,558,859,311,318,312,343,710,273,259,1274,625,245,235,250,NA,303,250,403,441,1507,388,530,107,868,NA,753,62,346,478,292) ylab = 'y' xlab = 'x' #'GNU S' R Code compiled by R2WASP v. 1.2.291 () #Author: root #To cite this work: Wessa P., (2012), Pearson Correlation (v1.0.6) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_correlation.wasp/ #Source of accompanying publication: Office for Research, Development, and Education # x <- x[!is.na(y)] y <- y[!is.na(y)] y <- y[!is.na(x)] x <- x[!is.na(x)] postscript(file="/var/wessaorg/rcomp/tmp/1iki51352119461.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) histx <- hist(x, plot=FALSE) histy <- hist(y, plot=FALSE) maxcounts <- max(c(histx$counts, histx$counts)) xrange <- c(min(x),max(x)) yrange <- c(min(y),max(y)) nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE) par(mar=c(4,4,1,1)) plot(x, y, xlim=xrange, ylim=yrange, xlab=xlab, ylab=ylab) par(mar=c(0,4,1,1)) barplot(histx$counts, axes=FALSE, ylim=c(0, maxcounts), space=0) par(mar=c(4,0,1,1)) barplot(histy$counts, axes=FALSE, xlim=c(0, maxcounts), space=0, horiz=TRUE) dev.off() lx = length(x) makebiased = (lx-1)/lx varx = var(x)*makebiased vary = var(y)*makebiased corxy <- cor.test(x,y,method='pearson', na.rm = T) cxy <- as.matrix(corxy$estimate)[1,1] #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/wessaorg/rcomp/createtable") a<-table.start() a<-table.row.start(a) a<-table.element(a,'Pearson Product Moment Correlation - Ungrouped Data',3,TRUE) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'Statistic',1,TRUE) a<-table.element(a,'Variable X',1,TRUE) a<-table.element(a,'Variable Y',1,TRUE) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/arithmetic_mean.htm','Mean',''),header=TRUE) a<-table.element(a,mean(x)) a<-table.element(a,mean(y)) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/biased.htm','Biased Variance',''),header=TRUE) a<-table.element(a,varx) a<-table.element(a,vary) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/biased1.htm','Biased Standard Deviation',''),header=TRUE) a<-table.element(a,sqrt(varx)) a<-table.element(a,sqrt(vary)) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/covariance.htm','Covariance',''),header=TRUE) a<-table.element(a,cov(x,y),2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/pearson_correlation.htm','Correlation',''),header=TRUE) a<-table.element(a,cxy,2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/coeff_of_determination.htm','Determination',''),header=TRUE) a<-table.element(a,cxy*cxy,2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,hyperlink('http://www.xycoon.com/ttest_statistic.htm','T-Test',''),header=TRUE) a<-table.element(a,as.matrix(corxy$statistic)[1,1],2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'p-value (2 sided)',header=TRUE) a<-table.element(a,(p2 <- as.matrix(corxy$p.value)[1,1]),2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'p-value (1 sided)',header=TRUE) a<-table.element(a,p2/2,2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'Degrees of Freedom',header=TRUE) a<-table.element(a,lx-2,2) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'Number of Observations',header=TRUE) a<-table.element(a,lx,2) a<-table.row.end(a) a<-table.end(a) table.save(a,file="/var/wessaorg/rcomp/tmp/2df1n1352119461.tab") try(system("convert tmp/1iki51352119461.ps tmp/1iki51352119461.png",intern=TRUE))