1 / 25

Spatial statistics in practice

Lab #3: - Jacobian estimation - estimating autoregressive models - estimating spatial filter models - estimating random effects models - mapping spatial filters. Spatial statistics in practice Center for Tropical Ecology and Biodiversity, Tunghai University & Fushan Botanical Garden.

cade
Download Presentation

Spatial statistics in practice

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lab #3:- Jacobian estimation- estimating autoregressive models- estimating spatial filter models- estimating random effects models - mapping spatial filters Spatial statistics in practice Center for Tropical Ecology and Biodiversity, Tunghai University & Fushan Botanical Garden

  2. SAS code for the Jacobian approximation FILENAME EIGEN 'D:\JYU-SUMMERSCHOOL2006\LAB#3\PR-EIG.TXT'; TITLE 'MATRIX C OR W JACOBIAN APPROXIMATION, IRREGULAR LATTICE'; *********************************************** * THE APPROPRIATE EIGENVALES MUST BE SELECTED * ***********************************************; DATA STEP0; INFILE EIGEN; INPUT ID LAMBDAC LAMBDAW; LAMBDA=LAMBDAW; X0=1; RUN; PROC UNIVARIATE NOPRINT; VAR LAMBDA; OUTPUT OUT=EXTREMEL MAX=LMAX MIN=LMIN; RUN; PROC PRINT DATA=EXTREMEL; VAR LMAX LMIN; RUN; PROC MEANS DATA=STEP0 NOPRINT; VAR X0; OUTPUT OUT=STEP0A SUM=N; RUN; DATA EXTREMEL(REPLACE=YES); SET EXTREMEL; SET STEP0A; FINISH = 0.999/LMAX; START = 0.999/LMIN; INC = (FINISH - START)/201; RUN; DATA STEP1; IF _N_=1 THEN SET EXTREMEL; SET STEP0; ARRAY JACOB{202} JAC1-JAC202; RHO = START; DO I = 1 TO 202; JACOB{I} = -LOG(1 - RHO*LAMBDA); RHO = RHO + INC; END; RUN; PROC MEANS NOPRINT; VAR JAC1-JAC202; OUTPUT OUT=JACOB1 MEAN=; RUN;

  3. PROC TRANSPOSE OUT=JACOB2; VAR JAC1-JAC202; RUN; DATA STEP2 (REPLACE=YES); SET EXTREMEL; DO RHO = START TO FINISH BY INC; OUTPUT; END; DROP START FINISH INC; RUN; DATA STEP2 (REPLACE=YES); SET JACOB2; SET STEP2; J = COL1; DROP COL1; RUN; PROC GPLOT; PLOT J*RHO; RUN; PROC NLIN DATA=STEP2 NOITPRINT MAXITER=15000 METHOD=MARQUARDT; PARMS A1=0.15 A2=0.15 D1=1.1 D2=1.1; BOUNDS D1<2 , D2<2 ; MODEL J = A1*(LOG(1 - RHO*LMIN)/(RHO*LMIN) + 1 - D1*LOG(1 - RHO*LMIN))+ A2*(LOG(1 - RHO*LMAX)/(RHO*LMAX) + 1 - D2*LOG(1 - RHO*LMAX)); OUTPUT OUT=TEMP3 ESS=ESS2 P=JHAT PARMS=A1 A2 D1 D2; RUN; PROC MEANS; VAR A1 A2 D1 D2; RUN; PROC GPLOT; PLOT J*RHO JHAT*RHO/OVERLAY; RUN; PROC UNIVARIATE DATA=STEP2 NOPRINT; VAR J; OUTPUT OUT=TOTALSS CSS=TSS; RUN; DATA TEMP3(REPLACE=YES); SET TOTALSS; SET TEMP3; RESS = ESS2/TSS; RUN; PROC PRINT; VAR ESS2 TSS RESS; RUN; one of three possible approximation equations

  4. SAS code for simultaneous autoregressive (SAR) modeling FILENAME ATTRIBUT 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-DEM&QUAD-DATA.TXT'; FILENAME CONN 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-CON.TXT'; FILENAME EIGEN 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-EIG.TXT'; TITLE 'SAR FOR PUERTO RICO DEM'; DATA STEP1; INFILE ATTRIBUT; INPUT IDDEM MELEV SELEV U V QUAD NAME$; U=U/1000; V=V/1000; Y=LOG(MELEV+17.5); * Y=(SELEV-25)**0.5; RUN; PROC SORT OUT=STEP1(REPLACE=YES); BY IDDEM; RUN; PROC STANDARD MEAN=0 STD=1 OUT=STEP1(REPLACE=YES); VAR Y; RUN; DATA STEP1 (REPLACE=YES); SET STEP1; INFILE CONN; INPUT MUNUM2 C1-C73; ARRAY CONY{73} CY1-CY73; ARRAY CON{73} C1-C73; CSUM = 0; DO I=1 TO 73; CSUM = CSUM + CON{I}; CONY{I} = Y*CON{I}; END; RUN; PROC MEANS SUM; VAR CSUM; RUN; PROC PRINT; VAR NAME IDDEM MUNUM2; RUN; PROC MEANS DATA=STEP1 NOPRINT; VAR CY1-CY73; OUTPUT OUT=CYOUT1 SUM=CY1-CY73; RUN; PROC TRANSPOSE DATA=CYOUT1 PREFIX=WY OUT=CYOUT2; VAR CY1-CY73; RUN; DATA STEP1 (REPLACE=YES); SET STEP1; SET CYOUT2; WY = WY1/CSUM; RUN; PROC REG; MODEL Y=WY; RUN;

  5. eigenvalues for the exact Jacobian DATA STEP2; INFILE EIGEN; INPUT IDE LAMBDAC LAMBDAW; LAMBDA=LAMBDAW; RUN; PROC TRANSPOSE DATA=STEP2 PREFIX=TLAM OUT=EOUT1; VAR LAMBDAW; RUN; DATA STEP1; SET STEP1; IF _N_ = 1 THEN SET EOUT1; RUN; PROC NLIN DATA=STEP1 NOITPRINT METHOD=MARQUARDT; PARMS RHO=0.5 B0=0; BOUNDS -1.5<RHO<1; ARRAY LAMBDAJ{73} TLAM1-TLAM73; JACOB = 0; DERJ = 0; DO I=1 TO 73; JACOB = JACOB + LOG(1 - RHO*LAMBDAJ{I}); DERJ = DERJ + -LAMBDAJ{I}/(1 - RHO*LAMBDAJ{I}); END; J=EXP(JACOB/73); DERJ = -DERJ/73; ZY = Y/J; MODEL ZY = (RHO*WY + B0*(1 - RHO))/J; OUTPUT OUT=TEMP1 PRED=YHAT R=YRESID PARMS=RHO B0; DER.RHO = ((RHO*WY + B0*(1 - RHO) - Y)*DERJ + WY - B0)/J; RUN; PROC REG; MODEL Y=YHAT; RUN; PROC UNIVARIATE NORMAL; VAR YRESID; RUN; SAS calculus not completely correct

  6. from the Jacobian approximation ************************** * * * JACOBIAN APPROXIMATION * * * **************************; PROC MEANS NOPRINT DATA=STEP2; VAR LAMBDA; OUTPUT OUT=EXTREMES MIN=LMIN MAX=LMAX; RUN; DATA STEP1(REPLACE=YES); SET STEP1; IF _N_=1 THEN SET EXTREMES(KEEP=LMIN LMAX); RUN; PROC NLIN DATA=STEP1 NOITPRINT METHOD=MARQUARDT MAXITER=500; PARMS RHO=0.5 B0=0; BOUNDS -1<RHO<1; A1 = 0.4276827; A2 = 0.2998183; D1 = 0.9853629; D2 = 1.0664149; IF RHO=0 THEN J=1; ELSE J = EXP(A1*(LOG(1 - RHO*LMIN)/(RHO*LMIN) + 1 - D1*LOG(1 - RHO*LMIN))+ A2*(LOG(1 - RHO*LMAX)/(RHO*LMAX) + 1 - D2*LOG(1 - RHO*LMAX)) ); IF RHO=0 THEN DERJ=0; ELSE DERJ = A1*(-LOG(1 - RHO*LMIN)/(LMIN*RHO**2) + (RHO*D1*LMIN - 1)/(RHO*(1 - RHO*LMIN)) ) + A2*(-LOG(1 - RHO*LMAX)/(LMAX*RHO**2) + (RHO*D2*LMAX - 1)/(RHO*(1 - RHO*LMAX)) ) ; ZY = Y*J; MODEL ZY = (RHO*WY + B0*(1 - RHO))*J; OUTPUT OUT=TEMP1 PRED=YHAT R=YRESID; DER.RHO = ( ( (RHO*WY + B0*(1 - RHO) ) - Y)*DERJ + WY - B0)*J; RUN;

  7. The AR-SAR model specification 4th order effect

  8. SAS code for spatial filter: Gaussian RV FILENAME ATTRIBUT 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-DEM&QUAD-DATA.TXT'; FILENAME MC 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-MC-EIG.TXT'; FILENAME CONN 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-CON.TXT'; FILENAME EVEC 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-EIGENVECTORS.TXT'; FILENAME OUTFILE 'D:\JYU-SUMMERSCHOOL2006\LAB#3\GAUSSIAN-SF.TXT'; OPTIONS LINESIZE=72; TITLE 'SPATIAL FILTER MODEL FOR PUERTO RICO: GAUSSIAN RV'; DATA STEP1; INFILE ATTRIBUT; INPUT IDDEM MELEV SELEV U V QUAD NAME$; * Y=LOG(MELEV+17.5); Y=(SELEV-25)**0.5; Y0=Y; RUN; PROC SORT DATA=STEP1 OUT=STEP1(REPLACE=YES); BY NAME; RUN; PROC UNIVARIATE NORMAL; VAR Y; RUN; PROC STANDARD MEAN=0 STD=1 OUT=STEP1(REPLACE=YES); VAR Y0; RUN; DATA STEP1 (REPLACE=YES); SET STEP1; INFILE EVEC LRECL=2048; INPUT IDE E1-E73; RUN; ***************************************************************** * * * THE SET OF CANDIDATE EIGENVECTORS IS DETERMINED BY MCADJ >= 0 * * * *****************************************************************; PROC REG OUTEST=COEF; MODEL Y = E1-E18/SELECTION=STEPWISE SLE=0.10; OUTPUT OUT=TEMP P=YHAT R=YRESID; RUN; PROC UNIVARIATE NORMAL; VAR YRESID; RUN; may wish to make a Bonferroni adjustment here stepwise regression

  9. PROC TRANSPOSE DATA=COEF PREFIX=B OUT=COEF2; VAR E1-E18; RUN; DATA COEF(REPLACE=YES); INFILE MC; INPUT ID LAM_MCM MC MCADJ; IF MCADJ<0.25 THEN DELETE; RUN; DATA COEF(REPLACE=YES); SET COEF; SET COEF2; IF B1='.' THEN DELETE; IF B1=0 THEN DELETE; BSQ=B1**2; EMC=BSQ*MC; P=1; RUN; PROC PRINT; RUN; PROC MEANS SUM NOPRINT; VAR LAM_MCM P EMC BSQ; OUTPUT OUT=EMC SUM=SUML P SPANUM SPADEN; RUN; PROC STANDARD DATA=TEMP MEAN=0 STD=1 OUT=TEMP(REPLACE=YES); VAR YRESID; RUN; DATA STEP2(REPLACE=YES); INFILE CONN LRECL=1024; INPUT ID C1-C73; RUN; DATA STEP2(REPLACE=YES); SET STEP2; SET TEMP(KEEP=YRESID); SET STEP1(KEEP=Y0); ARRAY CONN{73} C1-C73; ARRAY ZC{73} ZC1-ZC73; ARRAY Y0C{73} Y0C1-Y0C73; CSUM=0; DO I=1 TO 73; CSUM = CSUM + CONN{I}; ZC{I} = YRESID*CONN{I}; Y0C{I} = Y0*CONN{I}; END; Z=YRESID; X0=1; RUN; PROC MEANS SUM NOPRINT; VAR CSUM X0; OUTPUT OUT=CSUM SUM=CSUM N; RUN; PROC PRINT; VAR CSUM; RUN;

  10. PROC MEANS DATA=STEP2 NOPRINT; VAR Y0C1-Y0C73; OUTPUT OUT=Y0COUT1 SUM=Y0C1-Y0C73; RUN; PROC TRANSPOSE DATA=Y0COUT1 PREFIX=Y0C OUT=Y0COUT2; VAR Y0C1-Y0C73; RUN; PROC MEANS DATA=STEP2 NOPRINT; VAR ZC1-ZC73; OUTPUT OUT=ZCOUT1 SUM=ZC1-ZC73; RUN; PROC TRANSPOSE DATA=ZCOUT1 PREFIX=ZC OUT=ZCOUT2; VAR ZC1-ZC73; RUN; DATA STEP2(REPLACE=YES); SET STEP2(KEEP=X0 Z CSUM Y0); SET ZCOUT2(KEEP=ZC1); SET Y0COUT2(KEEP=Y0C1); Y0C=Y0C1; ZC=ZC1; DROP ZC1; RUN; PROC REG DATA=STEP2 OUTEST=DEN NOPRINT; MODEL CSUM=X0/NOINT; RUN; PROC REG DATA=STEP2 OUTEST=NUM0 NOPRINT; MODEL Y0C=Y0/NOINT; RUN; PROC REG DATA=STEP2 OUTEST=NUM NOPRINT; MODEL ZC=Z/NOINT; RUN; DATA STEP3; SET DEN; SET NUM0; SET NUM; SET CSUM(KEEP=CSUM N); SET EMC(KEEP=SUML SPANUM SPADEN P); MC0=Y0/X0; MC=Z/X0; EMC = -(1+(N/CSUM)*SUML)/(N-P-1); ZMC = (MC-EMC)/SQRT(2/CSUM); MCSPA=SPANUM/SPADEN; RUN; PROC PRINT; VAR MC0 MC EMC ZMC MCSPA; RUN; PROC STANDARD DATA=TEMP OUT=TEMP(REPLACE=YES) MEAN=0; VAR YHAT; RUN; DATA _NULL_; SET TEMP; FILE OUTFILE; PUT IDDEM Y YHAT; RUN; Residual diagnostics; MC distribution theory known for this case YHAT is the spatial filter; output for mapping purposes

  11. SAS code for spatial filter: binomial RV FILENAME INDATA 'D:\JYU-SUMMERSCHOOL2006\LAB#3\PR-POP-1899&2000.TXT'; FILENAME EVEC 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-EIGENVECTORS.TXT'; FILENAME OUTFILE 'D:\JYU-SUMMERSCHOOL2006\LAB#3\BINOMIAL-SF.TXT'; TITLE 'SPATIAL FILTER MODEL FOR PUERTO RICO: BINOMIAL RV'; DATA STEP1; INFILE INDATA; INPUT ID P1899 U1988 P2000 U2000 QUAD NAME$; Y=U2000/P2000; RUN; PROC SORT DATA=STEP1 OUT=STEP1(REPLACE=YES); BY NAME; RUN; DATA STEP1 (REPLACE=YES); SET STEP1; INFILE EVEC LRECL=2048; INPUT IDE E1-E73; RUN; ******************************************************************** * THE SET OF CANDIDATE EIGENVECTORS IS DETERMINED BY MCADJ >= 0.25 * ********************************************************************; PROC LOGISTIC; MODEL U2000/P2000 = E1-E18/SELECTION=STEPWISE SLE=0.10 SCALE=WILLIAMS; OUTPUT OUT=TEMP P=YHAT; RUN; PROC REG; MODEL Y=YHAT; RUN; PROC GENMOD; MODEL U2000/P2000 = E1 E3 E4 E10 E14 E18/DIST=BINOMIAL SCALE=DEVIANCE; OUTPUT OUT=SF XBETA=XBETA P=YHAT; RUN; DATA STEP1(REPLACE=YES); SET STEP1; Y=(U2000-1280)/(P2000-1057); TRY=LOG(Y/(1-Y)); W=1/((P2000-1057)*Y*(1-Y)); RUN; PROC REG; MODEL TRY=E1-E18/NOINT SELECTION=STEPWISE SLE=0.10; WEIGHT W; RUN; PROC STANDARD DATA=SF OUT=SF(REPLACE=YES) MEAN=0; VAR XBETA; RUN; DATA _NULL_; SET SF; FILE OUTFILE; YRESID=Y-YHAT; PUT ID Y XBETA YRESID; RUN; stepwise selection quasi- likelihood normal approximation XBETA is the spatial filter; output for mapping purposes

  12. SAS code for spatial filter: Poisson RV FILENAME INDATA1 'D:\JYU-SUMMERSCHOOL2006\LAB#3\PR-POP-1899&2000.TXT'; FILENAME INDATA2 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-AREAS-COMPETITION.TXT'; FILENAME EVEC 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-EIGENVECTORS.TXT'; FILENAME OUTFILE 'D:\JYU-SUMMERSCHOOL2006\LAB#3\POISSON-SF.TXT'; OPTIONS LINESIZE=72; TITLE 'SPATIAL FILTER MODEL FOR PUERTO RICO: POISSON RV'; DATA STEP1A; INFILE INDATA1; INPUT ID1 P1899 U1988 P2000 U2000 QUAD NAME$; DENOM=100000000000; RUN; PROC SORT DATA=STEP1A OUT=STEP1A(REPLACE=YES); BY NAME; RUN; DATA STEP1B; INFILE INDATA2; INPUT ID2 AREAM AREACT MCT_RATIO TRMCT_RATIO AREAPT MPT_RATIO TRMPT_RATIO NAME$; LNAREAM=LOG(AREAM); RUN; PROC SORT DATA=STEP1B OUT=STEP1B(REPLACE=YES); BY NAME; RUN; DATA STEP1; MERGE STEP1A STEP1B; BY NAME; RUN; DATA STEP1(REPLACE=YES); SET STEP1; INFILE EVEC LRECL=2048; INPUT IDE E1-E73; TRY=LOG(P2000/AREAM - 441830); RUN; ******************************************************************** * * * THE SET OF CANDIDATE EIGENVECTORS IS DETERMINED BY MCADJ >= 0.25 * * * ********************************************************************; PROC GENMOD; MODEL P2000 = E1-E18/DIST=POISSON OFFSET=LNAREAM; RUN; denominator for binomial approximation offset variable must be converted to its log form no stepwise Poisson regression option in SAS

  13. quasi- likelihood ***************************************************************************** * * * ORDER OF REMOVAL (BACKWARD ELIMINATION): 13, 5, 17, 14, 11, 15, 16, 9, 18 * * * *****************************************************************************; PROC GENMOD; MODEL P2000 = E1-E4 E6-E8 E10 E12/DIST=POISSON OFFSET=LNAREAM SCALE=DEVIANCE; RUN; ********************************************************************* * * * ORDER OF REMOVAL (BACKWARD ELIMINATION): 17, 14, 11, 13, 15, 5, 9 * * * *********************************************************************; PROC GENMOD; MODEL P2000 = E1-E4 E6-E8 E10 E12 E16 E18/DIST=NB OFFSET=LNAREAM; OUTPUT OUT=TEMP XBETA=XBETA P=YHAT; RUN; DATA TEMP(REPLACE=YES); SET TEMP; Y=P2000/AREAM; YHAT=YHAT/AREAM; RUN; PROC REG; MODEL Y=YHAT; RUN; PROC LOGISTIC; MODEL P2000/DENOM = E1-E18/SELECTION=BACKWARD OFFSET=LNAREAM SCALE=WILLIAMS SLS=0.10; RUN; PROC REG; MODEL TRY=E1-E18/SELECTION=BACKWARD SLS=0.10; RUN; PROC STANDARD DATA=TEMP OUT=TEMP(REPLACE=YES) MEAN=0; VAR XBETA; RUN; DATA _NULL_; SET TEMP; FILE OUTFILE; XBETA=XBETA-LNAREAM; YRESID=Y-YHAT; PUT ID1 Y XBETA YRESID; RUN; negative binomial binomial approximation normal approximation XBETA is the spatial filter; output for mapping purposes

  14. ArcView mapping of spatial filters Dark red: very high Light red: high Gray: medium Light green: low Dark green: very low • The SAS output files need to be converted to *.dbf files with column headers • MCs can be computed with MapStat3 for the binomial and Poisson regressions residuals Poisson binomial Gaussian

  15. SAS code for spatially structured random effects linear modeling FILENAME INDATA1 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-AREAS-COMPETITION.TXT'; FILENAME INDATA2 'D:\JYU-SUMMERSCHOOL2006\LAB#1\PR-DEM&QUAD-DATA.TXT'; TITLE 'MIXED MODEL FOR THE PR DATA'; ************************************************************************** * READ IN GEOREFERENCED DATA;THEN CENTER THE SELEDTED ATTRIBUTE VARIABLE * **************************************************************************; DATA STEP1; INFILE INDATA1 LRECL=1024; INPUT ID AREAM AREACT MCT_RATIO TRMCT_RATIO AREAPT MPT_RATIO TRMPT_RATIO NAME$; * Y = LOG(AREAM + 0.001); * Y = TRMPT_RATIO; RUN; PROC SORT DATA=STEP1 OUT=STEP1(REPLACE=YES); BY NAME; RUN; DATA STEP2; INFILE INDATA2 LRECL=1024; INPUT IDDEM MELEV SELEV U V QUAD NAME$; U=U/1000; V=V/1000; * Y=LOG(MELEV+17.5); Y=(SELEV-25)**0.5; X0=1; RUN; PROC SORT DATA=STEP2 OUT=STEP2(REPLACE=YES); BY NAME; RUN; DATA STEP3; MERGE STEP1 STEP2; BY NAME; RUN; PROC REG; MODEL Y=/; RUN; PROC MIXED METHOD=ML MAXITER=500 COVTEST SCORING=10; MODEL Y=X0/NOINT S; REPEATED/SUB=INTERCEPT; RUN;

  16. maximum likelihood estimation PROC MIXED METHOD=ML MAXITER=500 COVTEST SCORING=10; MODEL Y=X0/NOINT S; PARMS (13, 10); REPEATED/SUB=INTERCEPT TYPE=SP(EXP) (U V); RUN; PROC MIXED METHOD=ML MAXITER=500 COVTEST SCORING=10; MODEL Y=X0/NOINT S; PARMS (13, 10); REPEATED/SUB=INTERCEPT TYPE=SP(GAU) (U V); RUN; PROC MIXED METHOD=ML MAXITER=500 COVTEST SCORING=10; MODEL Y=X0/NOINT S; PARMS (13, 10); REPEATED/SUB=INTERCEPT TYPE=SP(SPH) (U V); RUN; PROC MIXED METHOD=ML MAXITER=500 COVTEST SCORING=10; MODEL Y=X0/NOINT S; PARMS (0.1, 13); REPEATED/SUB=INTERCEPT TYPE=SP(POW) (U V); RUN; PROC MIXED METHOD=ML MAXITER=500 COVTEST SCORING=10; MODEL Y=X0/NOINT S; PARMS (13, 10, 1); REPEATED/SUB=INTERCEPT TYPE=SP(MATERN) (U V); RUN; initial parameter values semivariogram model geocodings random effect

  17. OLS: 12.5690 Linear mixed model (Gaussian assumption) without nugget

  18. OLS: 12.5690 Linear mixed model (Gaussian assumption) with nugget

  19. SAS code for spatially structured random effects generalized linear modeling FILENAME INDATA 'D:\JYU-SUMMERSCHOOL2006\LAB#3\PR-URBAN-DATA.TXT'; FILENAME EVECS 'D:\JYU-SUMMERSCHOOL2006\LAB#3\PR-EXPANDED-EVECS.TXT'; OPTIONS LINESIZE=72; TITLE 'GENERALIZED LINEAR MIXED MODEL FOR THE PR DATA'; ******************************************************************* * READ IN GEOREFERENCED DATA; THEN CENTER THE ELEVATION VARIABLE * *******************************************************************; DATA STEP1; INFILE INDATA; INPUT ID Y1899 Y1910 Y1920 Y1930 Y1935 Y1940 Y1950 Y1960 Y1970 Y1980 Y1990 Y2000 MELEV SELEV U V AAR NAME$; NUM=Y1899; U=U/1000; V=V/1000; ELEV=LOG(MELEV/SELEV-0.6); UTELEV=MELEV/SELEV-0.6; IF AAR=1 THEN ISJ=1; ELSE ISJ=0; IF AAR=2 THEN IA =1; ELSE IA =0; IF AAR=3 THEN IM =1; ELSE IM =0; IF AAR=4 THEN IP =1; ELSE IP =0; IF AAR=5 THEN IC =1; ELSE IC =0; IA=IA-ISJ; IM=IM-ISJ; IP=IP-ISJ; IC=IC-ISJ; DENOM=100; Y=NUM/DENOM; P=Y; RUN;

  20. PROC STANDARD MEAN=0 STD=1 OUT=STEP1(REPLACE=YES); VAR ELEV; RUN; DATA STEP2; INFILE EVECS LRECL=1024; INPUT IDE E1-E45; RUN; DATA STEP2(REPLACE=YES); SET STEP2; SET STEP1; RUN; DATA STEP4; SET STEP2(KEEP=ID DENOM ELEV UTELEV IA IP IM IC Y1899 E1-E45); NUM=Y1899; TIME=1899; DROP Y1899; RUN; DATA STEP3; SET STEP2(KEEP=ID DENOM ELEV UTELEV IA IP IM IC Y1910 E1-E45); NUM=Y1910; TIME=1910; DROP Y1910; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1920 E1-E45); NUM=Y1920; TIME=1920; DROP Y1920; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1930 E1-E45); NUM=Y1930; TIME=1930; DROP Y1930; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1935 E1-E45); NUM=Y1935; TIME=1935; DROP Y1935; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1940 E1-E45); NUM=Y1940; TIME=1940; DROP Y1940; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1950 E1-E45); NUM=Y1950; TIME=1950; DROP Y1950; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1960 E1-E45); NUM=Y1960; TIME=1960; DROP Y1960; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1970 E1-E45); NUM=Y1970; TIME=1970; DROP Y1970; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1980 E1-E45); NUM=Y1980; TIME=1980; DROP Y1980; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y1990 E1-E45); NUM=Y1990; TIME=1990; DROP Y1990; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; DATA STEP3(REPLACE=YES); SET STEP2(KEEP=ID DENOM ELEV UTELEV IP IM IC Y2000 E1-E45); NUM=Y2000; TIME=2000; DROP Y2000; RUN; DATA STEP4(REPLACE=YES); SET STEP4 STEP3; RUN; PROC SORT OUT=STEP5; BY ID DESCENDING TIME; RUN;

  21. PROC NLMIXED DATA=STEP5; PARMS B0=0 BE=0 BIM=0 BIP=0 BIC=0 S2U=1 BE1 =0 BE2 =0 BE3 =0 BE4 =0 BE5 =0 BE6 =0 BE7 =0 BE9 =0 BE10=0 BE11=0 BE16=0 BE17=0 BE18=0 BE19=0 BE20=0 BE21=0 BE23=0 BE26=0 BE27=0 BE33=0 BE35=0; ETA = B0 + ALPHA + BE*ELEV + BIM*IM + BIP*IP + BIC*IC + BE1*E1+BE2*E2+BE3*E3+BE4*E4+BE5*E5+BE6*E6+BE7*E7+BE9*E9+BE10*E10+BE11*E11+BE16*E16+ BE17*E17+BE18*E18+BE19*E19+BE20*E20+BE21*E21+BE23*E23+BE26*E26+BE27*E27+BE33*E33+BE35*E35; P = EXP(ETA)/(1+EXP(ETA)); MODEL NUM ~ BINOMIAL(100,P); RANDOM ALPHA ~ NORMAL(0,S2U) SUBJECT=ID; PREDICT P OUT=POUT; RUN; PROC NLMIXED DATA=STEP5; PARMS B0=0 BE=0 BIM=0 BIC=0 S2U=1 BE2 =0 BE6 =0 BE10=0 BE11=0 BE19=0 BE35=0; ETA = B0 + ALPHA + BE*ELEV + BIM*IM + BIC*IC + BE2*E2+ BE6*E6+ BE10*E10+BE11*E11+ BE19*E19+ BE35*E35; P = EXP(ETA)/(1+EXP(ETA)); MODEL NUM ~ BINOMIAL(100,P); RANDOM ALPHA ~ NORMAL(0,S2U) SUBJECT=ID; PREDICT P OUT=POUT; RUN; PROC SORT OUT=POUT(REPLACE=YES); BY TIME ID; RUN; DATA POUT(REPLACE=YES); SET POUT; SET STEP4(KEEP=NUM); P=NUM/100; RUN; PROC REG; MODEL P=PRED; BY TIME; RUN; PROC GPLOT; PLOT P*PRED; BY TIME; RUN; PROC GENMOD DATA=STEP2; MODEL NUM/DENOM= ELEV E1 E4 E6 E7 E11 E17 E23 E26 E33/DIST=B SCALE=P; OUTPUT OUT=POUT P=PRED; RUN;

  22. Spatial filters structuring the random effects MC = 0.274 GR = 0.795 MC = 0.877 GR = 0.227 MC = -0.459 GR = 1.582

  23. The random effect MC = -0.074, GR = 0.973

  24. What you should have learned in today’s lab: • Calculate the J approximation for a given connectivity matrix (both C and W) • Normal approximation: estimate AR, SAR, AR-SAR and SF models • Estimate a SF GLM (Poisson or binomial) • Estimate a univariate LMM • Estimate a bivariate (time) SF-GLMM • Construct maps of the SFs • Construct maps of the random effects

More Related