\> \# File: CDF_CHEM.R

\> \# Purpose: Population estimation for Region 10 and states for
Chemistry

\> \# Data files include final design information and

\> \# only have one record for each site from the probability sample

\> \# Programmer: Tony Olsen

\> \# Date: October 11, 2005

\> \# Modified: Dan McKenzie

\> \# Date: January 27, 2006

\> \# Data required: chem_wgt3.csv

\> \#

\> \# Libraries required: psurvey.analysis 2.9

\> \# R scripts required:

\> \# cont.ash.r

\> \# ash1.wgt.r

\> \# cont.plot.r

\> \# plot.cont.r

\> \# plotcdf.r

\>

\> \# load psurvey.analysis

\> require(\'psurvey.analysis\')

\[1\] TRUE

\> \# source required R scripts

\> source(\'plot.cont.r\')

\> source(\'plotcdf.r\')

\> source(\'cont.plot.r\')

\> source(\'cont.ash.r\')

\> source(\'ash1.wgt.r\')

\>

\>

\> \#### Example Population estimates for Chemistry continuous variables

\> \# Read in status-weight-data file

\> chem \<- read.csv(\'Original Data/chem_wgt3.csv\')

\> head(chem)

SITE_ID VIS_STAT STATE EPAREG EVALUATE USED REPEAT LG_RIVER

1 WIDP99-0501 PRIMARY ID 10 Y Y

2 WIDP99-0502 PRIMARY ID 10 Y Y

3 WIDP99-0503 PRIMARY ID 10 Y Y

4 WIDP99-0504 PRIMARY ID 10 Y Y

5 WIDP99-0508 PRIMARY ID 10 Y Y

6 WIDP99-0510 PRIMARY ID 10 Y Y

RF3NAME RF3ID COUNTY MAP24K

1 17060303 603 0.00 IDAHO Cayuse Junction

2 GRIMES CR 17050112 41 12.97 BOISE Grimes Pass

3 17060308 984 0.00 LATAH McGary Butte

4 COEUR D\'ALENE R- N FK 17010301 48 18.75 KOOTENAI Spades Mountain

5 17050113 711 0.00 ELMORE Ross Peak

6 RATTLE CREEK 17010213 122 0.00 BONNER Trestle Peak

MAP100K MAP250K ECO ECOREGL3_DES ECOL3NAME_DES ECO10_DES ECO3_DES

1 MISSOULA WES HAMILTON 15 16 Idaho Batholith MT-NROCK MT

2 DEADWOOD RIV CHALLIS 15 16 Idaho Batholith MT-NROCK MT

3 POTLATCH PULLMAN 15 15 Northern Rockies MT-NROCK MT

4 COEUR D\'ALEN SPOKANE 15 15 Northern Rockies MT-NROCK MT

5 IDAHO CITY HAILEY 15 16 Idaho Batholith MT-NROCK MT

6 SANDPOINT SANDPOINT 15 15 Northern Rockies MT-NROCK MT

STRAHLER LAT_DD LONG_DD XLAT_DD XLON_DD POPLAT_DD POPLONG_DD OROCODE

1 2 46.56095 -114.7808 46.56306 -114.7806 46.56306 -114.7806 humid

2 2 44.03402 -115.8135 44.03278 -115.8192 44.03278 -115.8192 humid

3 2 46.75859 -116.3564 46.75861 -116.3569 46.75861 -116.3569 humid

4 3 47.79687 -116.5089 47.79750 -116.5097 47.79750 -116.5097 humid

5 1 43.68952 -115.0850 43.68944 -115.0858 43.68944 -115.0858 humid

6 2 48.32913 -116.1561 48.32861 -116.1564 48.32861 -116.1564 humid

PARTITON STRATUM WGT_CATY PANEL PANEL_ADJ OVERSAMP DIVISION STATEWID
CATY_N

1 1 17 162 0 0 0 2 1 4

2 1 17 162 0 0 0 2 1 4

3 1 17 162 0 0 0 2 1 4

4 1 17 163 0 0 0 2 1 3

5 2 17 161 0 0 0 2 1 4

6 2 17 162 0 0 0 2 1 4

INIT_WGT REMAP REMAP_N REMAP_WT SPECIAL SPEC_N SPEC_WGT SPECREMP SPR_N

1 3622.419 0 0 0 NA NA NA NA NA

2 3622.419 0 0 0 NA NA NA NA NA

3 3622.419 0 0 0 NA NA NA NA NA

4 2897.936 0 0 0 NA NA NA NA NA

5 7244.839 0 0 0 NA NA NA NA NA

6 3622.419 0 0 0 NA NA NA NA NA

SPR_WGT PLAIN PLAIN_N PLAIN_WT WGTCLASS YEAR VISIT_NO DATE_COL FLOWPAT

1 NA NA NA NA IDSW 2000 1 9/4/2000 P

2 NA NA NA NA IDSW 2000 1 8/23/2000 P

3 NA NA NA NA IDSW 2000 1 9/7/2000 P

4 NA NA NA NA IDSW 2000 1 9/12/2000 P

5 NA NA NA NA IDSW 2000 1 9/13/2000 P

6 NA NA NA NA IDSW 2000 1 9/10/2000 P

SITESAMP LOC_NAME COM_IM INDEXVIS XSTATUS VALXSTAT

1 Y WEST FORK PAPOOSE Y SAMPLEABLE WADEABLE

2 Y GRIMES CREEK Y SAMPLEABLE WADEABLE

3 Y NNT LONG MEADOWS CREEK Y SAMPLEABLE WADEABLE

4 Y LITTLE NORTH FORK COEUR D\'ALENE Y SAMPLEABLE WADEABLE

5 Y BADGER CREEK Y SAMPLEABLE WADEABLE

6 Y RATTLE CREEK Y SAMPLEABLE WADEABLE

TARGCLAS TARGCLAS2 STATUS_1 SURVEY ST_ORDER EC_ORD ECOREGL3

1 CANDIDATE TARGET-SAMPLED TS ST_WIDE 2nd HUMID 2nd 16

2 CANDIDATE TARGET-SAMPLED TS ST_WIDE 2nd HUMID 2nd 16

3 CANDIDATE TARGET-SAMPLED TS ST_WIDE 2nd HUMID 2nd 15

4 CANDIDATE TARGET-SAMPLED TS ST_WIDE 3rd HUMID 3rd 15

5 CANDIDATE TARGET-SAMPLED TS ST_WIDE 0-1st HUMID 1st 16

6 CANDIDATE TARGET-SAMPLED TS ST_WIDE 2nd HUMID 2nd 15

ECO10 ECO3 ORIG_WGT WGT_EXTENT WGT_COND SAMPCHEM PHSTVL NO3

1 MT-NROCK MT 3622.419 571.0931 716.2289 Yes 8.07 11.42

2 MT-NROCK MT 3622.419 571.0931 716.2289 Yes 7.74 0.70

3 MT-NROCK MT 3622.419 571.0931 716.2289 Yes 7.55 0.70

4 MT-NROCK MT 2897.936 456.8745 572.9831 Yes 7.50 2.14

5 MT-NROCK MT 7244.839 1142.1862 1432.4578 Yes 7.65 5.71

6 MT-NROCK MT 3622.419 571.0931 716.2289 Yes 7.03 24.27

\> dim(chem)

\[1\] 293 75

\> names(chem)

\[1\] \"SITE_ID\" \"VIS_STAT\" \"STATE\" \"EPAREG\"

\[5\] \"EVALUATE\" \"USED\" \"REPEAT\" \"LG_RIVER\"

\[9\] \"RF3NAME\" \"RF3ID\" \"COUNTY\" \"MAP24K\"

\[13\] \"MAP100K\" \"MAP250K\" \"ECO\" \"ECOREGL3_DES\"

\[17\] \"ECOL3NAME_DES\" \"ECO10_DES\" \"ECO3_DES\" \"STRAHLER\"

\[21\] \"LAT_DD\" \"LONG_DD\" \"XLAT_DD\" \"XLON_DD\"

\[25\] \"POPLAT_DD\" \"POPLONG_DD\" \"OROCODE\" \"PARTITON\"

\[29\] \"STRATUM\" \"WGT_CATY\" \"PANEL\" \"PANEL_ADJ\"

\[33\] \"OVERSAMP\" \"DIVISION\" \"STATEWID\" \"CATY_N\"

\[37\] \"INIT_WGT\" \"REMAP\" \"REMAP_N\" \"REMAP_WT\"

\[41\] \"SPECIAL\" \"SPEC_N\" \"SPEC_WGT\" \"SPECREMP\"

\[45\] \"SPR_N\" \"SPR_WGT\" \"PLAIN\" \"PLAIN_N\"

\[49\] \"PLAIN_WT\" \"WGTCLASS\" \"YEAR\" \"VISIT_NO\"

\[53\] \"DATE_COL\" \"FLOWPAT\" \"SITESAMP\" \"LOC_NAME\"

\[57\] \"COM_IM\" \"INDEXVIS\" \"XSTATUS\" \"VALXSTAT\"

\[61\] \"TARGCLAS\" \"TARGCLAS2\" \"STATUS_1\" \"SURVEY\"

\[65\] \"ST_ORDER\" \"EC_ORD\" \"ECOREGL3\" \"ECO10\"

\[69\] \"ECO3\" \"ORIG_WGT\" \"WGT_EXTENT\" \"WGT_COND\"

\[73\] \"SAMPCHEM\" \"PHSTVL\" \"NO3\"

\>

\> \# add equal area coordinates for local variance estimation

\> tmp \<- marinus(chem\$POPLAT_DD, chem\$POPLONG_DD)

\> chem\$xmarinus \<- tmp\[,\'x\'\]

\> chem\$ymarinus \<- tmp\[,\'y\'\]

\>

\>

\> \# Set up data frames for cdf estimation

\> sites.chem \<- data.frame(SiteID=chem\$SITE_ID,

\+ Use=(chem\$STATUS_1==\'TS\' & chem\$SAMPCHEM==\'Yes\'))

\> subpop.chem \<- data.frame(SiteID=chem\$SITE_ID,

\+ Region10=rep(\'Region 10\', nrow(chem)),

\+ State=chem\$STATE

\+ )

\> dsgn.chem \<- data.frame(siteID=chem\$SITE_ID,

\+ stratum=chem\$STRATUM,

\+ wgt=chem\$WGT_COND,

\+ xcoord=chem\$xmarinus,

\+ ycoord=chem\$ymarinus )

\>

\> data.cont.chem \<- data.frame(siteID=chem\$SITE_ID,

\+ pH=chem\$PHSTVL,

\+ Nitrate=chem\$NO3

\+ )

\>

\>

\> \# Do cdf on sites

\> chemcdf \<- cont.analysis(sites=sites.chem,

\+ subpop=subpop.chem,

\+ design=dsgn.chem,

\+ data.cont=data.cont.chem,

\+ pctval=c(5, 10, 25, 50, 75, 90, 95)

\+ )

\>

\> write.table(chemcdf\$Pct,\'chem_pct.csv\',sep = \",\",col.names=NA)

\>

\>

\> \## Estimate empirical density function for chemistry variables

\> \# define ind.type frame for use in both density and plot functions

\> ind.type.chem \<- rep(\'Continuous\',2)

\>

\> \# below is for all variables

\> ind.type.chem \<- rep(\'Continuous\',20)

\>

\> \# define ab matrix for use in density function

\> \# note: set lower limit=0 for any variables that you will

\> \# be plotting using a logrithmic scale

\> ab.chem \<- matrix(c(3,14,

\+ 0,700),

\+ nc=2,

\+ byrow=TRUE)

\>

\>

\> \# ash density estimates

\> density.chem \<- cont.ash(sites=sites.chem,

\+ subpop=subpop.chem,

\+ design=dsgn.chem,

\+ data.cont=data.cont.chem,

\+ nbin=100,

\+ ab=ab.chem,

\+ ind.type=ind.type.chem,

\+ vartype = \"Local\",

\+ )

\>

\> #define x-axis labels

\> xlabels.chem \<- c(\'pH\', \'Nitrate (ueq/L)\')

\> names(xlabels.chem) \<-c(\'pH\',\'Nitrate\')

\>

\>

\>

\>

\> \# plot all subpopulations and indicators

\> cont.plot(pdffile=\'chemistry_plots.pdf\',chemcdf,
denlist=density.chem,

\+ ind.type=ind.type.chem, figlab=\'chem\',

\+ xlab=xlabels.chem, ylab.u=\'Stream Length (km)\',

\+ logx=c(FALSE, TRUE),

\+ width=7, height=8)

\>

\>

\>