R version 3.2.2 (2015-08-14) -- "Fire Safety" Copyright (C) 2015 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(1747,1245,1182,958,1000,1044,875,939,736,905,796,372,1326,668,962,912,1119,891,931,1047,982,1098,714,128,1784,828,1199,1095,977,1338,975,840,1324,1236,883,177,2186,809,1434,1365,1247,1476,1211,990,1205,1238,952,204,2135,1157,1290,1071,1169,1431,945,1034,1100,1297,921,236,1990,966,1326,908,1206,1861,929,1296,1332,1352,1040,148,2090,1435,1124,1319,1436,1774,1566,1385,1147,1274,625,52,1990,1154,954,887,825,966,954,770,1838,1371,589,116,1898,712,1175,1240,1329,1550,1201,938,1030,1060,1035,635,2565,910,1304,1331,1681,1983,1021,1061,1292,1274,1024,568,2570,1125,1600,1492,2492,3523,990,869,1310,979,1244,442,2956,1055,2004,1462,1144,1454,1538,1388,1547,1570,1535,1352,1888,999,1158,1342,1443,1519,1267,1454,987,1430,1254,734) > 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] 156 > (np <- floor(n / par1)) [1] 13 > 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] [,11] [,12] [,13] [1,] 1747 1326 1784 2186 2135 1990 2090 1990 1898 2565 2570 2956 1888 [2,] 1245 668 828 809 1157 966 1435 1154 712 910 1125 1055 999 [3,] 1182 962 1199 1434 1290 1326 1124 954 1175 1304 1600 2004 1158 [4,] 958 912 1095 1365 1071 908 1319 887 1240 1331 1492 1462 1342 [5,] 1000 1119 977 1247 1169 1206 1436 825 1329 1681 2492 1144 1443 [6,] 1044 891 1338 1476 1431 1861 1774 966 1550 1983 3523 1454 1519 [7,] 875 931 975 1211 945 929 1566 954 1201 1021 990 1538 1267 [8,] 939 1047 840 990 1034 1296 1385 770 938 1061 869 1388 1454 [9,] 736 982 1324 1205 1100 1332 1147 1838 1030 1292 1310 1547 987 [10,] 905 1098 1236 1238 1297 1352 1274 1371 1060 1274 979 1570 1430 [11,] 796 714 883 952 921 1040 625 589 1035 1024 1244 1535 1254 [12,] 372 128 177 204 236 148 52 116 635 568 442 1352 734 > 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] 983.2500 898.1667 1054.6667 1193.0833 1148.8333 1196.1667 1268.9167 [8] 1034.5000 1150.2500 1334.5000 1553.0000 1583.7500 1289.5833 > arr.sd [1] 328.2956 299.4491 387.2623 465.5475 430.9448 474.0796 524.0086 511.8129 [9] 343.8562 530.7260 878.4892 491.9655 299.6251 > arr.range [1] 1375 1198 1607 1982 1899 1842 2038 1874 1263 1997 3081 1901 1154 > (lm1 <- lm(arr.sd~arr.mean)) Call: lm(formula = arr.sd ~ arr.mean) Coefficients: (Intercept) arr.mean -141.7413 0.4977 > (lnlm1 <- lm(log(arr.sd)~log(arr.mean))) Call: lm(formula = log(arr.sd) ~ log(arr.mean)) Coefficients: (Intercept) log(arr.mean) -2.235 1.175 > (lm2 <- lm(arr.range~arr.mean)) Call: lm(formula = arr.range ~ arr.mean) Coefficients: (Intercept) arr.mean -137.896 1.594 > postscript(file="/var/wessaorg/rcomp/tmp/11nwg1448230617.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/2s2p61448230617.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/387jn1448230617.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/4vx051448230617.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/5nth31448230617.tab") > > try(system("convert tmp/11nwg1448230617.ps tmp/11nwg1448230617.png",intern=TRUE)) character(0) > try(system("convert tmp/2s2p61448230617.ps tmp/2s2p61448230617.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 0.822 0.192 1.016