LAB 2: EXERCISE 1:Direct Product Operator EXERCISE 2: J Function EXERCISE 3:Sample variogram in 4 directions EXERCISE 4:Semivariogram of residuals (after fitting mean with proc reg) EXERCISE 5:Theoretical and sample semivariogram (fitting by eye) EXERCISE 6:Semivariogram WITH CASTNET DATA /* FYI: ABOUT THE HELP COMMAND: help: e.g. we want to learn about proc print: the HELP PRINT needs to be in the toolbox or you can use the pull down menu HELP> SAS SYSTEM> then for print REPORT WRITING > PRINT or you can use the PGM editor if you enter the statement dm 'help print'; run; and then submit this also brings up the SAS HELP */ EXERCISE 1 ********** /* Direct Product Operator: @ takes the direct product of two matrices matrix1@matrix2 The direct product operator (@) produces a new matrix that is the direct product (also called the Kronecker product) of matrix1 and matrix2, usually denoted by .The number of rows in the new matrix equals the product of the number of rows in matrix1 and the number of rows in matrix2; the number of columns in the new matrix equals the product of the number of columns in matrix1 and the number of columns in matrix2. For example, the statements */ proc iml; a={1 2, 3 4}; b={0 2}; c=a@b; print c; run; /* result in C 2 rows 4 cols (numeric) 0 2 0 4 0 6 0 8 */ **The statement proc iml; a={1 2, 3 4}; b={0 2}; d=b@a; print d; run; /*results in D 2 rows 4 cols (numeric) 0 0 2 4 0 0 6 8 */ ********************************************************************** EXERCISE 2 ********** /* J Function creates a matrix of identical values J( nrow<, ncol<, value>>) The inputs to the J function are as follows: nrow is a numeric matrix or literal giving the number of rows. ncol is a numeric matrix or literal giving the number of columns. value is a numeric or character matrix or literal for filling the rows and columns of the matrix. The J function creates a matrix with nrow rows and ncol columns with all elements equal to value. If ncol is not specified, it defaults to nrow. If value is not specified, it defaults to 1. The REPEAT and SHAPE functions can also perform this operation, and they are more general. Examples of the J function are shown below:*/ proc iml; b=j(3); print b; run; /* B 3 rows 3 cols (numeric) 1 1 1 1 1 1 1 1 1 */ proc iml; r=j(5,2,'xyz'); print r; run; /* R 5 rows 2 cols (character, size 3) xyz xyz xyz xyz xyz xyz xyz xyz xyz xyz */ **HOW TO GET A MATRIX FROM THE DATASET soilph.dat (using the j function) data phdat; infile '/ncsu/st733_info/www/soilph.dat'; input row col value; proc iml; use phdat; /* defines current data to read from*/ read all var {row col value}; nval=nrow(row); /* number of values in matrix*/ nrow=max(row); /* number of row */ ncol=max(col); /* number of col */ matval=j(nrow,ncol,0); /* initialize matrix */ /* put values in corresponding position in matrix matval */ do i=1 to nval; matval[row[i,1],col[i,1]]=value[i,1]; end; print matval;run; ************************************************************************ EXERCISE 3 *********** *** Sample variogram in 4 directions. options nodate nocenter nonumber ps=67 ls=120 ; data soilph; infile '/ncsu/st733_info/www/soilph.dat'; input row col ph; proc variogram data=soilph outv=outv; compute lagd=0.7 maxlag=8 ndir=4; /* compute lagd=0.7 maxlag=8 ndir=4 robust; (to get robust semivgram.)*/ coordinates xc=col yc=row; run; proc print data=outv; run; proc sort data=outv; by angle distance; goptions horigin=1in; proc gplot data=outv; where (distance > .); plot variog *distance=angle / vaxis =axis2 haxis=axis1 nolegend ; symbol1 c=black i=join l=1 h=3 v='0'; symbol2 c=red i=join l=2 h=3 v='1'; symbol3 c=green i=join l=3 h=3 v='2'; symbol4 c=blue i=join l=4 h=3 v='3'; axis1 minor=none label=(c=black 'lag distance') offset=(3,3); axis2 minor=(number =1) label=(c=black 'Vgram' ) offset=(3,3); title 'Sample variogram in 4 directions'; footnote '0--0 0 deg, 1--1 45 deg, 2--2 90 deg, 3--3 135 deg.'; run; ***************************************************************** Output of exercise 3: OBS VARNAME ANGLE ATOL LAG COUNT DISTANCE AVERAGE VARIOG COVAR 1 PH . . -1 121 . 4.48463 . 0.045759 2 PH 0 22.5 0 0 . . . . 3 PH 0 22.5 1 110 1.00000 4.49555 0.028760 0.017286 4 PH 0 22.5 2 0 . . . . 5 PH 0 22.5 3 99 2.00000 4.50379 0.036326 0.009589 6 PH 0 22.5 4 88 3.00000 4.51236 0.051354 -0.004628 7 PH 0 22.5 5 160 3.16228 4.51966 0.053866 -0.004692 8 PH 0 22.5 6 217 4.07942 4.50878 0.058865 -0.009702 9 PH 0 22.5 7 186 5.06388 4.50134 0.068428 -0.018524 10 PH 0 22.5 8 108 5.38516 4.50144 0.068762 -0.019062 11 PH 45 22.5 0 0 . . . . 12 PH 45 22.5 1 0 . . . . 13 PH 45 22.5 2 100 1.41421 4.50173 0.038371 0.010411 14 PH 45 22.5 3 180 2.23607 4.50511 0.046708 0.001360 15 PH 45 22.5 4 81 2.82843 4.50580 0.049707 -0.002209 16 PH 45 22.5 5 144 3.60555 4.50151 0.055896 -0.009352 17 PH 45 22.5 6 190 4.39483 4.49572 0.061169 -0.013851 18 PH 45 22.5 7 112 5.00000 4.48993 0.059218 -0.011977 19 PH 45 22.5 8 145 5.77212 4.48022 0.064001 -0.017295 20 PH 90 22.5 0 0 . . . . 21 PH 90 22.5 1 110 1.00000 4.48932 0.029931 0.017344 22 PH 90 22.5 2 0 . . . . 23 PH 90 22.5 3 99 2.00000 4.48652 0.035297 0.010937 24 PH 90 22.5 4 88 3.00000 4.47736 0.030749 0.013158 25 PH 90 22.5 5 160 3.16228 4.48975 0.036201 0.008000 26 PH 90 22.5 6 217 4.07942 4.48879 0.039159 0.006041 27 PH 90 22.5 7 186 5.06388 4.49073 0.041852 0.002406 28 PH 90 22.5 8 108 5.38516 4.50264 0.041273 0.003241 29 PH 135 22.5 0 0 . . . . 30 PH 135 22.5 1 0 . . . . 31 PH 135 22.5 2 100 1.41421 4.50075 0.032662 0.015400 32 PH 135 22.5 3 180 2.23607 4.50382 0.032087 0.015062 33 PH 135 22.5 4 81 2.82843 4.51185 0.030292 0.016778 34 PH 135 22.5 5 144 3.60555 4.51472 0.037060 0.009146 35 PH 135 22.5 6 190 4.39483 4.51667 0.040318 0.006155 36 PH 135 22.5 7 112 5.00000 4.52116 0.040981 0.004847 37 PH 135 22.5 8 145 5.77212 4.52495 0.046633 -0.000395 ********************************************************************** EXERCISE 4 *********** **Semivariogram of residuals: options nodate nocenter nonumber ps=67 ls=120 ; data soilph; infile '/ncsu/st733_info/www/soilph.dat'; input row col ph; r2=row*row; c2=col*col; rc=row*col; run; proc reg data=soilph; model ph=row col r2 rc c2; output out=resdat r=resid; *proc print data=resdat; run; proc variogram data=resdat outv=outv; compute lagd=0.7 maxlag=8 ndir=4; coordinates xc=col yc=row; var resid; run; proc print data=outv; run; proc sort data=outv; by angle distance; goptions horigin=1in; proc gplot data=outv; where (distance > .); plot variog *distance=angle / vaxis =axis2 haxis=axis1 nolegend ; symbol1 c=black i=join l=1 h=3 v='0'; symbol2 c=red i=join l=2 h=3 v='1'; symbol3 c=green i=join l=3 h=3 v='2'; symbol4 c=blue i=join l=4 h=3 v='3'; axis1 minor=none label=(c=black 'lag distance') offset=(3,3); axis2 minor=(number =1) label=(c=black 'Vgram' ) offset=(3,3); title 'Sample variogram in 4 directions'; footnote '0--0 0 deg, 1--1 45 deg, 2--2 90 deg, 3--3 135 deg.'; run; ************************************************************************* Output of exercise 4: Sample variogram in 4 directions Model: MODEL1 Dependent Variable: PH Analysis of Variance Sum of Mean Source DF Squares Square F Value Prob>F Model 5 1.27645 0.25529 6.966 0.0001 Error 115 4.21466 0.03665 C Total 120 5.49111 Root MSE 0.19144 R-square 0.2325 Dep Mean 4.48463 Adj R-sq 0.1991 C.V. 4.26880 Parameter Estimates Parameter Standard T for H0: Variable DF Estimate Error Parameter=0 Prob > |T| INTERCEP 1 4.197752 0.10797728 38.876 0.0001 ROW 1 0.119825 0.02642918 4.534 0.0001 COL 1 0.026320 0.02642918 0.996 0.3214 R2 1 -0.008831 0.00197057 -4.481 0.0001 RC 1 -0.004640 0.00174036 -2.666 0.0088 C2 1 -0.000363 0.00197057 -0.184 0.8540 Sample variogram in 4 directions OBS VARNAME ANGLE ATOL LAG COUNT DISTANCE AVERAGE VARIOG COVAR 1 RESID . . -1 121 . 0.000000 . 0.035122 2 RESID 0 22.5 0 0 . . . . 3 RESID 0 22.5 1 110 1.00000 -0.002329 0.028079 0.006558 4 RESID 0 22.5 2 0 . . . . 5 RESID 0 22.5 3 99 2.00000 -0.001445 0.032681 0.002110 6 RESID 0 22.5 4 88 3.00000 0.005653 0.039826 -0.004215 7 RESID 0 22.5 5 160 3.16228 0.012406 0.040832 -0.004169 8 RESID 0 22.5 6 217 4.07942 0.006137 0.041853 -0.004885 9 RESID 0 22.5 7 186 5.06388 0.009005 0.042745 -0.005361 10 RESID 0 22.5 8 108 5.38516 0.008600 0.041849 -0.005525 11 RESID 45 22.5 0 0 . . . . 12 RESID 45 22.5 1 0 . . . . 13 RESID 45 22.5 2 100 1.41421 0.004466 0.036761 -0.000408 14 RESID 45 22.5 3 180 2.23607 0.005181 0.041587 -0.005540 15 RESID 45 22.5 4 81 2.82843 0.004362 0.039756 -0.004312 16 RESID 45 22.5 5 144 3.60555 0.001624 0.040792 -0.005688 17 RESID 45 22.5 6 190 4.39483 -0.000186 0.041992 -0.005936 18 RESID 45 22.5 7 112 5.00000 -0.001461 0.034635 0.001808 19 RESID 45 22.5 8 145 5.77212 -0.002969 0.036666 -0.000344 20 RESID 90 22.5 0 0 . . . . 21 RESID 90 22.5 1 110 1.00000 0.004145 0.029906 0.005812 22 RESID 90 22.5 2 0 . . . . 23 RESID 90 22.5 3 99 2.00000 0.001039 0.034700 -0.000431 24 RESID 90 22.5 4 88 3.00000 -0.008179 0.029706 0.003102 25 RESID 90 22.5 5 160 3.16228 -0.009033 0.033462 -0.000577 26 RESID 90 22.5 6 217 4.07942 -0.005111 0.035589 -0.001529 27 RESID 90 22.5 7 186 5.06388 -0.002751 0.036770 -0.002905 28 RESID 90 22.5 8 108 5.38516 -0.002897 0.033850 0.000167 29 RESID 135 22.5 0 0 . . . . 30 RESID 135 22.5 1 0 . . . . 31 RESID 135 22.5 2 100 1.41421 0.001170 0.032400 0.003023 32 RESID 135 22.5 3 180 2.23607 -0.000751 0.030906 0.003523 33 RESID 135 22.5 4 81 2.82843 0.001130 0.029108 0.005159 34 RESID 135 22.5 5 144 3.60555 0.000914 0.032524 0.001075 35 RESID 135 22.5 6 190 4.39483 0.001418 0.033570 0.000227 36 RESID 135 22.5 7 112 5.00000 0.001924 0.032977 0.000292 37 RESID 135 22.5 8 145 5.77212 0.006167 0.035324 -0.001877 ******************************************************************* EXERCISE 5 *********** **Theoretical and sample semivariogram. proc variogram data=resdat outv=outo; compute lagd=0.7 maxlag=8 robust; coordinates xc=col yc=row; var resid; run; data outo2; set outo; th0= 0.008; th1= 0.03; th2=1/.8; vari= th0+ th1*(1-exp(-distance*th2) ); type='Exponential'; output; vari= variog; type='regular'; output; *proc print data=outo2; run; proc gplot data=outo2; where (distance > .); plot vari*distance=type; symbol1 c=black i=join l=1 v='0'; symbol2 c=red i=join l=2 v='1'; axis1 minor=none label=(c=black 'lag distance') offset=(3,3); axis2 minor=(number=1) label=(c=green 'semi-variogram') offset=(3,3); title 'Sample variogram versus theoretical'; footnote ' '; run; /* data outo3; set outo; /*Spherical variogram fit*/ th0=0.003; th1= 0.035; th2=2; if ( distance < 2) then vari= th0+ th1*( 3*distance/(2*th2)- ((distance/th2)**3)/2); if ( distance > 2) then vari= th1+th0; type='Spherical'; output; vari= variog; type='sample'; output; run; proc print data=outo3; run; data outo4; set outo; /*Gaussian variogram fit*/ th0=0.003; th1= 0.035; th2=1/2; vari= th0+ th1*( 1- exp(-distance*distance*th2*th2)); type='Gaussian'; output; vari= variog; type='sample'; output; run; proc print data=outo4; run; */ ************************************** Output from spherical fit: Sample variogram versus theoretical OBS VARNAME LAG COUNT DISTANCE AVERAGE VARIOG COVAR RVARIO TH0 TH1 TH2 VARI TYPE 1 RESID -1 121 . .0000000 . 0.035122 . .003 0.035 2 . Spherical 2 RESID -1 121 . .0000000 . 0.035122 . .003 0.035 2 . sample 3 RESID 0 0 . . . . . .003 0.035 2 . Spherical 4 RESID 0 0 . . . . . .003 0.035 2 . sample 5 RESID 1 220 1.00000 .0009081 0.028993 0.006185 0.027615 .003 0.035 2 0.027063 Spherical 6 RESID 1 220 1.00000 .0009081 0.028993 0.006185 0.027615 .003 0.035 2 0.028993 sample 7 RESID 2 200 1.41421 .0028181 0.034580 0.001307 0.041132 .003 0.035 2 0.033936 Spherical 8 RESID 2 200 1.41421 .0028181 0.034580 0.001307 0.041132 .003 0.035 2 0.034580 sample 9 RESID 3 558 2.15230 .0013569 0.035340 -0.000353 0.034009 .003 0.035 2 0.038000 Spherical 10 RESID 3 558 2.15230 .0013569 0.035340 -0.000353 0.034009 .003 0.035 2 0.035340 sample 11 RESID 4 338 2.91777 .0006585 0.034606 -0.000087 0.034516 .003 0.035 2 0.038000 Spherical 12 RESID 4 338 2.91777 .0006585 0.034606 -0.000087 0.034516 .003 0.035 2 0.034606 sample 13 RESID 5 608 3.37225 .0014887 0.036915 -0.002341 0.039947 .003 0.035 2 0.038000 Spherical 14 RESID 5 608 3.37225 .0014887 0.036915 -0.002341 0.039947 .003 0.035 2 0.036915 sample 15 RESID 6 814 4.22667 .0005614 0.038282 -0.003042 0.039504 .003 0.035 2 0.038000 Spherical 16 RESID 6 814 4.22667 .0005614 0.038282 -0.003042 0.039504 .003 0.035 2 0.038282 sample 17 RESID 7 596 5.03987 .0020389 0.037521 -0.002185 0.042475 .003 0.035 2 0.038000 Spherical 18 RESID 7 596 5.03987 .0020389 0.037521 -0.002185 0.042475 .003 0.035 2 0.037521 sample 19 RESID 8 506 5.60694 .0021335 0.036787 -0.001780 0.036119 .003 0.035 2 0.038000 Spherical 20 RESID 8 506 5.60694 .0021335 0.036787 -0.001780 0.036119 .003 0.035 2 0.036787 sample ************************************************************************ EXERCISE 6 ********** *WITH CASTNET DATA (LON, LAT, OZONE) options nodate nocenter nonumber ps=67 ls=120 ; data soilph; infile '/ncsu/st733_info/www/castnet.oz'; input lon lat oz; proc variogram data=soilph outv=outv; compute lagd=1.5 maxlag=8 ndir=4; /* compute lagd=1.5 maxlag=8 ndir=4 robust; (to get robust semivgram.)*/ coordinates xc=lon yc=lat; run; proc print data=outv; run; proc sort data=outv; by angle distance; goptions horigin=1in; proc gplot data=outv; where (distance > .); plot variog *distance=angle / vaxis =axis2 haxis=axis1 nolegend ; symbol1 c=black i=join l=1 h=3 v='0'; symbol2 c=red i=join l=2 h=3 v='1'; symbol3 c=green i=join l=3 h=3 v='2'; symbol4 c=blue i=join l=4 h=3 v='3'; axis1 minor=none label=(c=black 'lag distance') offset=(3,3); axis2 minor=(number =1) label=(c=black 'Vgram' ) offset=(3,3); title 'Sample variogram in 4 directions'; footnote '0--0 0 deg, 1--1 45 deg, 2--2 90 deg, 3--3 135 deg.'; run;