R version 3.3.1 (2016-06-21) -- "Bug in Your Hair" Copyright (C) 2016 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(3441,3431,3420,3400,3606,3596,3441,3338,3348,3348,3358,3379,3441,3472,3524,3565,3751,3730,3575,3348,3389,3431,3420,3472,3431,3503,3534,3544,3772,3730,3575,3348,3389,3358,3410,3524,3513,3493,3544,3575,3751,3761,3575,3307,3286,3348,3296,3462,3462,3400,3482,3534,3710,3761,3544,3286,3286,3203,3141,3276,3224,3100,3183,3255,3472,3555,3327,3162,3162,3100,3059,3141,3038,3017,3069,3141,3358,3400,3131,2935,2842,2749,2697,2800,2738,2749,2800,2842,3048,3079,2749,2594,2439,2335,2263,2366,2315,2397,2428,2459,2594,2676,2263,2160,1901,1736,1684,1829,1746,1850,1850,1860,1984,2067,1664,1519,1281,1126,1044,1250) > 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 Plot (v1.0.2) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_sdplot.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+1)) > ari <- array(0,dim=par1) > j <- 0 > for (i in 1:n) + { + j = j + 1 + ari[j] = ari[j] + 1 + arr[j,ari[j]] <- x[i] + if (j == par1) j = 0 + } > ari [1] 10 10 10 10 10 10 10 10 10 10 10 10 > arr [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 3441 3441 3431 3513 3462 3224 3038 2738 2315 1746 NA [2,] 3431 3472 3503 3493 3400 3100 3017 2749 2397 1850 NA [3,] 3420 3524 3534 3544 3482 3183 3069 2800 2428 1850 NA [4,] 3400 3565 3544 3575 3534 3255 3141 2842 2459 1860 NA [5,] 3606 3751 3772 3751 3710 3472 3358 3048 2594 1984 NA [6,] 3596 3730 3730 3761 3761 3555 3400 3079 2676 2067 NA [7,] 3441 3575 3575 3575 3544 3327 3131 2749 2263 1664 NA [8,] 3338 3348 3348 3307 3286 3162 2935 2594 2160 1519 NA [9,] 3348 3389 3389 3286 3286 3162 2842 2439 1901 1281 NA [10,] 3348 3431 3358 3348 3203 3100 2749 2335 1736 1126 NA [11,] 3358 3420 3410 3296 3141 3059 2697 2263 1684 1044 NA [12,] 3379 3472 3524 3462 3276 3141 2800 2366 1829 1250 NA > arr.sd <- array(NA,dim=par1) > arr.range <- array(NA,dim=par1) > arr.iqr <- array(NA,dim=par1) > for (j in 1:par1) + { + arr.sd[j] <- sqrt(var(arr[j,],na.rm=TRUE)) + arr.range[j] <- max(arr[j,],na.rm=TRUE) - min(arr[j,],na.rm=TRUE) + arr.iqr[j] <- quantile(arr[j,],0.75,na.rm=TRUE) - quantile(arr[j,],0.25,na.rm=TRUE) + } > overall.sd <- sqrt(var(x)) > overall.range <- max(x) - min(x) > overall.iqr <- quantile(x,0.75) - quantile(x,0.25) > postscript(file="/var/wessaorg/rcomp/tmp/1rszg1471293799.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(arr.sd,type='b',ylab='S.D.',main='Standard Deviation Plot',xlab='Periodic Index') > mtext(paste('# blocks = ',np)) > abline(overall.sd,0) > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/2prh61471293799.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(arr.range,type='b',ylab='range',main='Range Plot',xlab='Periodic Index') > mtext(paste('# blocks = ',np)) > abline(overall.range,0) > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/3yog41471293799.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(arr.iqr,type='b',ylab='IQR',main='Interquartile Range Plot',xlab='Periodic Index') > mtext(paste('# blocks = ',np)) > abline(overall.iqr,0) > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/4hy5q1471293799.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > z <- data.frame(t(arr)) > names(z) <- c(1:par1) > (boxplot(z,notch=TRUE,col='grey',xlab='Periodic Index',ylab='Value',main='Notched Box Plots - Periodic Subseries')) $stats [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] 1746.0 1850 1850.0 1860.0 2594 2676.0 1664 1519 1281 1126.0 1044 1250.0 [2,] 2738.0 2749 2800.0 2842.0 3048 3079.0 2749 2594 2439 2335.0 2263 2366.0 [3,] 3327.5 3250 3301.5 3327.5 3539 3575.5 3384 3224 3224 3151.5 3100 3208.5 [4,] 3441.0 3472 3524.0 3544.0 3751 3730.0 3575 3338 3348 3348.0 3358 3462.0 [5,] 3513.0 3503 3544.0 3575.0 3772 3761.0 3575 3348 3389 3431.0 3420 3524.0 $n [1] 10 10 10 10 10 10 10 10 10 10 10 10 $conf [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 2976.253 2888.76 2939.761 2976.753 3187.753 3250.234 2971.297 2852.268 [2,] 3678.747 3611.24 3663.239 3678.247 3890.247 3900.766 3796.703 3595.732 [,9] [,10] [,11] [,12] [1,] 2769.827 2645.365 2552.894 2660.895 [2,] 3678.173 3657.635 3647.106 3756.105 $out [1] 1984 2067 $group [1] 5 6 $names [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" Warning message: In bxp(list(stats = c(1746, 2738, 3327.5, 3441, 3513, 1850, 2749, : some notches went outside hinges ('box'): maybe set notch=FALSE > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/5hci51471293799.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > z <- data.frame(arr) > names(z) <- c(1:np) > (boxplot(z,notch=TRUE,col='grey',xlab='Block Index',ylab='Value',main='Notched Box Plots - Sequential Blocks')) $stats [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 3338 3348.0 3348.0 3286.0 3141 3059.0 2697.0 2263.0 1684.0 1044.0 NA [2,] 3353 3425.5 3399.5 3327.5 3281 3120.5 2821.0 2402.5 1865.0 1265.5 NA [3,] 3410 3472.0 3513.5 3503.0 3431 3172.5 3027.5 2743.5 2289.0 1705.0 NA [4,] 3441 3570.0 3559.5 3575.0 3539 3291.0 3136.0 2821.0 2443.5 1855.0 NA [5,] 3441 3751.0 3772.0 3761.0 3761 3472.0 3400.0 3079.0 2676.0 2067.0 NA $n [1] 12 12 12 12 12 12 12 12 12 12 0 $conf [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 3369.863 3406.093 3440.523 3390.114 3313.324 3094.734 2883.826 2552.619 [2,] 3450.137 3537.907 3586.477 3615.886 3548.676 3250.266 3171.174 2934.381 [,9] [,10] [,11] [1,] 2025.142 1436.125 NA [2,] 2552.858 1973.875 NA $out [1] 3606 3596 3555 $group [1] 1 1 6 $names [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" NA Warning message: In bxp(list(stats = c(3338, 3353, 3410, 3441, 3441, 3348, 3425.5, : some notches went outside hinges ('box'): maybe set notch=FALSE > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/6gqv71471293799.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > z <- data.frame(cbind(arr.sd,arr.range,arr.iqr)) > names(z) <- list('S.D.','Range','IQR') > (boxplot(z,notch=TRUE,col='grey',ylab='Overall Variability',main='Notched Box Plots')) $stats [,1] [,2] [,3] [1,] 558.4729 1653.0 570.750 [2,] 571.4386 1704.5 626.375 [3,] 611.9409 1808.5 648.625 [4,] 760.9065 2191.0 851.125 [5,] 822.3120 2376.0 971.000 $n [1] 12 12 12 $conf [,1] [,2] [,3] [1,] 525.5234 1586.604 546.115 [2,] 698.3585 2030.396 751.135 $out numeric(0) $group numeric(0) $names [1] "S.D." "Range" "IQR" Warning message: In bxp(list(stats = c(558.472878084116, 571.438566498107, 611.940947969509, : some notches went outside hinges ('box'): maybe set notch=FALSE > dev.off() null device 1 > > try(system("convert tmp/1rszg1471293799.ps tmp/1rszg1471293799.png",intern=TRUE)) character(0) > try(system("convert tmp/2prh61471293799.ps tmp/2prh61471293799.png",intern=TRUE)) character(0) > try(system("convert tmp/3yog41471293799.ps tmp/3yog41471293799.png",intern=TRUE)) character(0) > try(system("convert tmp/4hy5q1471293799.ps tmp/4hy5q1471293799.png",intern=TRUE)) character(0) > try(system("convert tmp/5hci51471293799.ps tmp/5hci51471293799.png",intern=TRUE)) character(0) > try(system("convert tmp/6gqv71471293799.ps tmp/6gqv71471293799.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.265 0.169 2.481