R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > x <- c(36439,36368,36290,36147,37615,37543,36439,35705,35777,35777,35848,35998,35998,35335,35043,35335,36368,36218,34822,33640,33419,32977,33276,33640,33497,33198,32614,33198,33718,33568,31873,31139,30405,29814,29743,30184,29593,29372,29151,30405,30548,29814,27826,26943,25547,24955,25247,25689,25689,25326,25247,26430,27385,26943,25468,24735,23189,22234,22968,23702,23702,22747,22676,23922,24735,24442,22968,22013,19947,19142,19434,20688,20759,18921,19584,21201,21935,21493,19506,18109,16492,15238,15751,16855,16563,14946,15459,17076,17960,17447,15459,14576,13251,11854,12075,13179,13322,11997,12218,14063,14504,13764,11042,9646,7801,5963,6554,7359,7217,5813,6625,8613,9496,9055,7288,5892,4417,2721,3021,3534) > par1 = '12' > par1 <- '12' > #'GNU S' R Code compiled by R2WASP v. 1.2.291 () > #Author: root > #To cite this work: Wessa P. (2012), Standard Deviation-Mean Plot (v1.0.6) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_smp.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > # > par1 <- as.numeric(par1) > (n <- length(x)) [1] 120 > (np <- floor(n / par1)) [1] 10 > arr <- array(NA,dim=c(par1,np)) > j <- 0 > k <- 1 > for (i in 1:(np*par1)) + { + j = j + 1 + arr[j,k] <- x[i] + if (j == par1) { + j = 0 + k=k+1 + } + } > arr [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 36439 35998 33497 29593 25689 23702 20759 16563 13322 7217 [2,] 36368 35335 33198 29372 25326 22747 18921 14946 11997 5813 [3,] 36290 35043 32614 29151 25247 22676 19584 15459 12218 6625 [4,] 36147 35335 33198 30405 26430 23922 21201 17076 14063 8613 [5,] 37615 36368 33718 30548 27385 24735 21935 17960 14504 9496 [6,] 37543 36218 33568 29814 26943 24442 21493 17447 13764 9055 [7,] 36439 34822 31873 27826 25468 22968 19506 15459 11042 7288 [8,] 35705 33640 31139 26943 24735 22013 18109 14576 9646 5892 [9,] 35777 33419 30405 25547 23189 19947 16492 13251 7801 4417 [10,] 35777 32977 29814 24955 22234 19142 15238 11854 5963 2721 [11,] 35848 33276 29743 25247 22968 19434 15751 12075 6554 3021 [12,] 35998 33640 30184 25689 23702 20688 16855 13179 7359 3534 > arr.mean <- array(NA,dim=np) > arr.sd <- array(NA,dim=np) > arr.range <- array(NA,dim=np) > for (j in 1:np) + { + arr.mean[j] <- mean(arr[,j],na.rm=TRUE) + arr.sd[j] <- sd(arr[,j],na.rm=TRUE) + arr.range[j] <- max(arr[,j],na.rm=TRUE) - min(arr[,j],na.rm=TRUE) + } > arr.mean [1] 36328.83 34672.58 31912.58 27924.17 24943.00 22201.33 18820.33 14987.08 [9] 10686.08 6141.00 > arr.sd [1] 641.5651 1227.9318 1574.7680 2144.7449 1627.7409 1958.4063 2323.3468 [8] 2065.5659 3110.8347 2334.2286 > arr.range [1] 1910 3391 3975 5593 5151 5593 6697 6106 8541 6775 > (lm1 <- lm(arr.sd~arr.mean)) Call: lm(formula = arr.sd ~ arr.mean) Coefficients: (Intercept) arr.mean 3171.66880 -0.05558 > (lnlm1 <- lm(log(arr.sd)~log(arr.mean))) Call: lm(formula = log(arr.sd) ~ log(arr.mean)) Coefficients: (Intercept) log(arr.mean) 12.7257 -0.5294 > (lm2 <- lm(arr.range~arr.mean)) Call: lm(formula = arr.range ~ arr.mean) Coefficients: (Intercept) arr.mean 9153.6126 -0.1654 > postscript(file="/var/wessaorg/rcomp/tmp/1c8j11344596654.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(arr.mean,arr.sd,main='Standard Deviation-Mean Plot',xlab='mean',ylab='standard deviation') > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/2kkad1344596654.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(arr.mean,arr.range,main='Range-Mean Plot',xlab='mean',ylab='range') > dev.off() null device 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,'Standard Deviation-Mean Plot',4,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Section',header=TRUE) > a<-table.element(a,'Mean',header=TRUE) > a<-table.element(a,'Standard Deviation',header=TRUE) > a<-table.element(a,'Range',header=TRUE) > a<-table.row.end(a) > for (j in 1:np) { + a<-table.row.start(a) + a<-table.element(a,j,header=TRUE) + a<-table.element(a,arr.mean[j]) + a<-table.element(a,arr.sd[j] ) + a<-table.element(a,arr.range[j] ) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/wessaorg/rcomp/tmp/3d2fu1344596654.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Regression: S.E.(k) = alpha + beta * Mean(k)',2,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'alpha',header=TRUE) > a<-table.element(a,lm1$coefficients[[1]]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'beta',header=TRUE) > a<-table.element(a,lm1$coefficients[[2]]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'S.D.',header=TRUE) > a<-table.element(a,summary(lm1)$coefficients[2,2]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'T-STAT',header=TRUE) > a<-table.element(a,summary(lm1)$coefficients[2,3]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'p-value',header=TRUE) > a<-table.element(a,summary(lm1)$coefficients[2,4]) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/wessaorg/rcomp/tmp/4wzh41344596654.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Regression: ln S.E.(k) = alpha + beta * ln Mean(k)',2,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'alpha',header=TRUE) > a<-table.element(a,lnlm1$coefficients[[1]]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'beta',header=TRUE) > a<-table.element(a,lnlm1$coefficients[[2]]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'S.D.',header=TRUE) > a<-table.element(a,summary(lnlm1)$coefficients[2,2]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'T-STAT',header=TRUE) > a<-table.element(a,summary(lnlm1)$coefficients[2,3]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'p-value',header=TRUE) > a<-table.element(a,summary(lnlm1)$coefficients[2,4]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Lambda',header=TRUE) > a<-table.element(a,1-lnlm1$coefficients[[2]]) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/wessaorg/rcomp/tmp/5fixn1344596654.tab") > > try(system("convert tmp/1c8j11344596654.ps tmp/1c8j11344596654.png",intern=TRUE)) character(0) > try(system("convert tmp/2kkad1344596654.ps tmp/2kkad1344596654.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.992 0.248 1.284