R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-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(152512,151944,151368,150176,161968,161344,152512,146640,147208,147208,147840,148976,150744,150744,149608,146640,161968,164304,160776,152512,156048,150744,153136,154280,155472,152512,153136,148976,161968,166072,162544,156048,163112,155472,162544,161968,163736,157240,164304,163736,174336,171944,162544,157808,164304,155472,161968,163112,165504,160208,163112,164880,171376,166072,159008,151368,158440,139000,148408,153704,159008,151368,151368,151368,155472,149608,141912,135472,140144,121904,133080,139576,140768,134272,134840,133080,139000,134840,126640,120712,130736,108968,123104,129544,129544,121904,114840,114272,120712,114840,103672,95976,104240,84808,102472,111872,114840,108344,100136,106008,108344,106576,88904,80704,86568,68904,87144,93640,98936,90104,81840,86568,88904,84232,66568,58872,65936,46504,67704,80704) > 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,] 152512 150744 155472 163736 165504 159008 140768 129544 114840 98936 [2,] 151944 150744 152512 157240 160208 151368 134272 121904 108344 90104 [3,] 151368 149608 153136 164304 163112 151368 134840 114840 100136 81840 [4,] 150176 146640 148976 163736 164880 151368 133080 114272 106008 86568 [5,] 161968 161968 161968 174336 171376 155472 139000 120712 108344 88904 [6,] 161344 164304 166072 171944 166072 149608 134840 114840 106576 84232 [7,] 152512 160776 162544 162544 159008 141912 126640 103672 88904 66568 [8,] 146640 152512 156048 157808 151368 135472 120712 95976 80704 58872 [9,] 147208 156048 163112 164304 158440 140144 130736 104240 86568 65936 [10,] 147208 150744 155472 155472 139000 121904 108968 84808 68904 46504 [11,] 147840 153136 162544 161968 148408 133080 123104 102472 87144 67704 [12,] 148976 154280 161968 163112 153704 139576 129544 111872 93640 80704 > 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] 151641.33 154292.00 158318.67 163375.33 158423.33 144190.00 129708.67 [8] 109929.33 95842.67 76406.00 > arr.sd [1] 5138.870 5452.686 5348.391 5472.115 9011.263 10743.845 8835.491 [8] 12314.689 13751.002 15218.567 > arr.range [1] 15328 17664 17096 18864 32376 37104 31800 44736 45936 52432 > (lm1 <- lm(arr.sd~arr.mean)) Call: lm(formula = arr.sd ~ arr.mean) Coefficients: (Intercept) arr.mean 24315.1031 -0.1132 > (lnlm1 <- lm(log(arr.sd)~log(arr.mean))) Call: lm(formula = log(arr.sd) ~ log(arr.mean)) Coefficients: (Intercept) log(arr.mean) 25.365 -1.386 > (lm2 <- lm(arr.range~arr.mean)) Call: lm(formula = arr.range ~ arr.mean) Coefficients: (Intercept) arr.mean 84786.7220 -0.3983 > postscript(file="/var/wessaorg/rcomp/tmp/1l4po1438333783.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/2etun1438333783.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/3lqje1438333783.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/4n8271438333783.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/5rw281438333783.tab") > > try(system("convert tmp/1l4po1438333783.ps tmp/1l4po1438333783.png",intern=TRUE)) character(0) > try(system("convert tmp/2etun1438333783.ps tmp/2etun1438333783.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.919 0.165 1.074