# File: ASSESSMENT_CLASS.R
# Purpose: condition class estimation for Region 10 and states
#          Useful for assessment document.
#          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: condclass_wgt3.csv

# Libraries required: psurvey.analysis 2.9

# Read in condition class data
assess <- read.csv('Original Data/condclass_wgt3.csv')
names(assess)

# Change Condition Class variables to factors to work with cat.analysis program
#   and treat missing value as missing, not as a class level
assess$BUG_COND <- as.factor(assess$BUG_COND)
assess$PTL_COND <- as.factor(assess$PTL_COND)
assess$NTL_COND <- as.factor(assess$NTL_COND)


########

## Need equal area coordinates for variance estimation
#  (uses x-site coords when available, design coords otherwise)
tmp <- marinus(assess$POPLAT_DD, assess$POPLONG_DD)
assess$xmarinus <- tmp[,'x']
assess$ymarinus <- tmp[,'y']

# Set up data for estimation
# which sites to use?
# create site file for biological/chem/habitat, invasive plant
# estimates

sites.cond <- data.frame(SiteID=assess$SITE_ID,
                     Use=(assess$USED=='Y'& !is.na(assess$WGT_COND)))

# want estimates for what subpopulations?
# Create subpop variables
subpop.cond <- data.frame(SiteID=assess$SITE_ID, 
                          Region10=rep('Region 10', nrow(assess)),
                          States=assess$STATE)

# Provide design information
#  -- change wgt variable for biology/habitat, 


dsgn.cond <- data.frame(siteID=assess$SITE_ID, 
                        stratum=assess$STRATUM, wgt=assess$WGT_COND, 
                        xcoord=assess$xmarinus, ycoord=assess$ymarinus)

# Provide categorical indicator data
data.cat.cond <- data.frame(siteID=assess$SITE_ID,
							 assess$BUG_COND,
							 assess$PTL_COND,
							 assess$NTL_COND)


# Do status population estimation for biology, nutrient, and habitat indicators
popstatus.cond <- cat.analysis(sites = sites.cond,
					subpop = subpop.cond, 
					design = dsgn.cond,
  				data.cat = data.cat.cond, 
					vartype = "Local", 
					conf = 95)

# write results out
write.table(popstatus.cond, 
         'condition_assessments.csv',
          sep = ",", col.names=NA)