Nothing
inst/doc/mastifVignette.R
In mastif: Mast Inference and Forecasting
## ----outform, echo=F----------------------------------------------------------
insertPlot <- function( file, caption ){
# outputFormat = knitr::opts_knit$get( "rmarkdown.pandoc.to" )
# if( outputFormat == 'latex' )
# paste( "![ ", caption, " ]( ", file, " )", sep="" )
}
bigskip <- function( ){
# outputFormat = knitr::opts_knit$get( "rmarkdown.pandoc.to" )
# if( outputFormat == 'latex' )
# "\\bigskip"
# else
"<br>"
}
## ----getFiles, eval = F, echo=F-----------------------------------------------
# Rcpp::sourceCpp( '../RcppFunctions/cppFns.cpp' )
# source( '../RFunctions/mastifFunctions.r' )
# library( RANN )
# library( robustbase )
## ----simSetup0, eval = F------------------------------------------------------
# seedNames <- specNames <- 'acerRubr'
# sim <- list( nplot=5, nyr=10, ntree=30, ntrap=40, specNames = specNames, seedNames = seedNames )
## ----sim0, eval = F-----------------------------------------------------------
# inputs <- mastSim( sim ) # simulate dispersal data
# seedData <- inputs$seedData # year, plot, trap, seed counts
# treeData <- inputs$treeData # year, plot, tree data
# xytree <- inputs$xytree # tree locations
# xytrap <- inputs$xytrap # trap locations
# formulaFec <- inputs$formulaFec # fecundity model
# formulaRep <- inputs$formulaRep # maturation model
# trueValues <- inputs$trueValues # true states and parameter values
## ----formulaFec, eval = F-----------------------------------------------------
# formulaFec
## ----treeData0, eval = F------------------------------------------------------
# head( treeData )
## ----xytree, eval = F---------------------------------------------------------
# head( xytree, 5 )
## ----treeData1, eval = F------------------------------------------------------
# head( seedData )
## ----map1a, eval = F----------------------------------------------------------
# dataTab <- table( treeData$plot, treeData$year )
#
# w <- which( dataTab > 0, arr.ind=T ) # a plot-year with observations
# w <- w[sample( nrow( w ), 1 ), ]
#
# mapYears <- as.numeric( colnames( dataTab )[w[2]] )
# mapPlot <- rownames( dataTab )[w[1]]
# inputs$mapPlot <- mapPlot
# inputs$mapYears <- mapYears
# inputs$treeSymbol <- treeData$diam
# inputs$SCALEBAR <- T
#
# mastMap( inputs )
## ----map3, eval=F-------------------------------------------------------------
# inputs$treeSymbol <- trueValues$fec
# inputs$treeScale <- 2
# inputs$trapScale <- 1
# mastMap( inputs )
## ----map4, eval=F-------------------------------------------------------------
# inputs$treeSymbol <- trueValues$repr
# inputs$treeScale <- .5
# mastMap( inputs )
## ----hist0, eval = F----------------------------------------------------------
# par( mfrow=c( 2, 1 ), bty='n', mar=c( 4, 4, 1, 1 ) )
# seedData <- inputs$seedData
# seedNames <- inputs$seedNames
#
# hist( as.matrix( seedData[, seedNames] ) , nclass=100,
# xlab = 'seed count', ylab='per trap', main='' )
# hist( trueValues$fec, nclass=100, xlab = 'seeds produced', ylab = 'per tree', main = '' )
## ----mast2, eval = F----------------------------------------------------------
# output <- mastif( inputs = inputs, ng = 4000, burnin = 500 )
## ----tabPars0, eval = F-------------------------------------------------------
# summary( output )
## ----pars, eval = F-----------------------------------------------------------
# plotPars <- list( trueValues = trueValues )
# mastPlot( output, plotPars )
## ---- eval=F------------------------------------------------------------------
# output <- mastif( inputs = output, ng = 5000, burnin = 3000 )
## ----restart, eval=F----------------------------------------------------------
# inputs$predList <- list( mapMeters = 3, plots = mapPlot, years = mapYears )
# output <- mastif( inputs = output, ng = 2000, burnin = 1000 )
## ----mapout, eval = F---------------------------------------------------------
# output$mapPlot <- mapPlot
# output$mapYears <- mapYears
# output$treeScale <- 1.5
# output$trapScale <- .8
# output$PREDICT <- T
# output$scaleValue <- 20
# output$plotScale <- 1
# output$COLORSCALE <- T
# output$LEGEND <- T
#
# mastMap( output )
## ----to rmd, eval=F-----------------------------------------------------------
# plotPars$RMD <- 'pdf'
# mastPlot( output, plotPars )
## ----simSetup, eval = F-------------------------------------------------------
# specNames <- c( 'pinuTaeda', 'pinuEchi', 'pinuVirg' )
# seedNames <- c( 'pinuTaeda', 'pinuEchi', 'pinuVirg', 'pinuUNKN' )
# sim <- list( nyr=4, ntree=25, nplot=10, ntrap=50, specNames = specNames,
# seedNames = seedNames )
## ----sim, eval = F------------------------------------------------------------
# inputs <- mastSim( sim ) # simulate dispersal data
# R <- inputs$trueValues$R # species to seedNames probability matrix
# round( R, 2 )
## ----mast3, eval = F----------------------------------------------------------
# output <- mastif( inputs = inputs, ng = 2000, burnin = 1000 )
## ----tabPars, eval = F--------------------------------------------------------
# summary( output )
## ----pars0, eval = F----------------------------------------------------------
# plotPars <- list( trueValues = inputs$trueValues )
# mastPlot( output, plotPars )
## ----again, eval=F------------------------------------------------------------
# tab <- with( inputs$seedData, table( plot, year ) )
# years <- as.numeric( colnames( tab )[tab[1, ] > 0] ) # years for 1st plot
# output$predList <- list( plots = 'p1', years = years )
# output <- mastif( inputs = output, ng = 3000, burnin = 1500 )
# mastPlot( output )
## ---- eval=F------------------------------------------------------------------
# specNames <- c( 'pinuTaed', 'pinuEchi' )
#
# #seeds never differentiated:
# seedNames <- c( 'pinuUNKN' )
#
# #one species sometimes differentiated:
# seedNames <- c( 'pinuTaed', 'pinuUNKN' )
#
# #both species sometimes differentiated:
# seedNames <- c( 'pinuTaed', 'pinuEchi', 'pinuUNKN' )
## ----map21, eval=F------------------------------------------------------------
# library( repmis )
# d <- "https://github.com/jimclarkatduke/mast/blob/master/liriodendronExample.rData?raw=True"
# repmis::source_data( d )
# mapList <- list( treeData = treeData, seedData = seedData,
# specNames = specNames, seedNames = seedNames,
# xytree = xytree, xytrap = xytrap, mapPlot = 'DUKE_BW',
# mapYears = 2011:2014, treeSymbol = treeData$diam,
# treeScale = .7, trapScale = 1.5, plotScale = 2,
# SCALEBAR=T, scaleValue=50 )
# mastMap( mapList )
## ----litu1, eval=F------------------------------------------------------------
# head( treeData, 3 )
## ----litu2, eval=F------------------------------------------------------------
# head( seedData, 3 )
## ----fit, eval=F--------------------------------------------------------------
# formulaFec <- as.formula( ~ diam ) # fecundity model
# formulaRep <- as.formula( ~ diam ) # maturation model
#
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData, xytree = xytree,
# xytrap = xytrap )
# output <- mastif( inputs, formulaFec, formulaRep, ng = 3000, burnin = 1000 )
## ----more, eval=F-------------------------------------------------------------
# output$predList <- list( mapMeters = 10, plots = 'DUKE_EW', years = 2010:2015 )
# output <- mastif( inputs = output, ng = 4000, burnin = 1000 )
## ----plotmydata1, eval=F------------------------------------------------------
# mastPlot( output )
## ----outpars, eval=F----------------------------------------------------------
# summary( output )
## ----fitSum, eval=F-----------------------------------------------------------
# output$fit
## ----nopred, eval=F-----------------------------------------------------------
# #group plots in regions for year effects
# region <- rep( 'sApps', nrow( treeData ) )
# region[ as.character( treeData$plot ) == 'DUKE_EW' ] <- 'piedmont'
#
# treeData$region <- region
#
# formulaFec <- as.formula( ~ diam )
# formulaRep <- as.formula( ~ diam )
# yearEffect <- list( groups = 'region' )
# randomEffect <- list( randGroups = 'treeID', formulaRan = as.formula( ~ 1 ) )
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData,
# xytree = xytree, xytrap = xytrap,
# priorDist = 28, priorVDist = 15, maxDist = 50, minDist = 15,
# minDiam = 25, maxF = 1e+6,
# randomEffect = randomEffect, yearEffect = yearEffect )
# output <- mastif( inputs, formulaFec, formulaRep, ng = 2000, burnin = 1000 )
# mastPlot( output )
## ---- yr, eval=F--------------------------------------------------------------
# yearEffect <- list( groups = c( 'species', 'region' ) ) # year effects
## ---- eval=F------------------------------------------------------------------
# inputs$priorList <- list( minDiam = 15, maxDiam = 60 )
## ---- eval=F------------------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/priorParametersPinus.csv?raw=True"
# download.file( d, destfile="priorParametersPinus.csv" )
#
# specNames <- c( "pinuEchi", "pinuRigi", "pinuStro", "pinuTaed", "pinuVirg" )
# seedNames <- c( "pinuEchi", "pinuRigi", "pinuStro", "pinuTaed", "pinuVirg", "pinuUNKN" )
# priorTable <- mastPriors( "priorParametersPinus.csv", specNames,
# code = 'code4', genus = 'pinus' )
## ---- eval=F, echo=F----------------------------------------------------------
#
# # THIS BLOCK IS OMITTED
#
# load( '../combinedSites/pinus.Rdata' )
#
# pl <- c( 'CWT_118', 'CWT_218', 'DUKE_BW', 'DUKE_EW', 'MARS_F', 'MARS_P' )
# treeData <- treeData[treeData$plot %in% pl, ]
# seedData <- seedData[seedData$plot %in% pl, ]
# xytree <- xytree[xytree$plot %in% pl, ]
# xytrap <- xytrap[xytrap$plot %in% pl, ]
#
# count <- as.matrix( seedData[, -c( 1:5 )] )
# count <- count[, colSums( count, na.rm=T ) > 0]
# count <- count[, -1]
# seedData <- cbind( seedData[, 1:5], count )
#
# treeData <- treeData[, 1:6]
#
# seedNames <- colnames( count )
# specNames <- sort( unique( treeData$species ) )
#
# save( treeData, seedData, xytree, xytrap,
# specNames, seedNames, file='pinusExample.rData' )
#
# #load( 'pinusExample.rData' )
## ----priorBeta, eval=F--------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/pinusExample.rdata?raw=True"
# repmis::source_data( d )
## ---- eval=F------------------------------------------------------------------
# formulaRep <- as.formula( ~ diam )
# formulaFec <- as.formula( ~ diam + I( diam^2 ) )
# betaPrior <- list( pos = 'diam', neg = 'I( diam^2 )' )
#
# inputs <- list( treeData = treeData, seedData = seedData, xytree = xytree,
# xytrap = xytrap, specNames = specNames, seedNames = seedNames,
# betaPrior = betaPrior, priorTable = priorTable,
# seedTraits = seedTraits )
# output <- mastif( inputs = inputs, formulaFec, formulaRep,
# ng = 500, burnin = 200 )
#
# # restart
# output <- mastif( output, formulaFec, formulaRep, ng = 5000, burnin = 1000 )
# mastPlot( output )
## ----fill, eval=F-------------------------------------------------------------
# # randomly remove years for this example:
# years <- sort( unique( treeData$year ) )
# sy <- sample( years, 5 )
# treeData <- treeData[treeData$year %in% sy, ]
# treeData[1:10, ]
## ---- removed, eval=F---------------------------------------------------------
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData,
# xytree = xytree, xytrap = xytrap, priorTable = priorTable,
# seedTraits = seedTraits )
# inputs <- mastFillCensus( inputs, beforeFirst=10, afterLast=10 )
# inputs$treeData[1:10, ]
## ---- treeYr table, eval=F----------------------------------------------------
# # original data
# table( treeData$year )
#
# # filled census
# table( inputs$treeData$year )
# table( seedData$year )
## ---- treeOnly, eval=F--------------------------------------------------------
# p <- 3
# inputs <- mastFillCensus( inputs, p = p )
# treeData <- inputs$treeData
## ----arr, eval=F--------------------------------------------------------------
# region <- c( 'CWT', 'DUKE', 'HARV' )
# treeData$region <- 'SCBI'
# for( j in 1:length( region ) ){
# wj <- which( startsWith( treeData$plot, region[j] ) )
# treeData$region[wj] <- region[j]
# }
# inputs$treeData <- treeData
# yearEffect <- list( groups = c( 'species', 'region' ), p = p )
## ----def, eval=F--------------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/def.csv?raw=True"
# download.file( d, destfile="def.csv" )
## ----types, eval=F------------------------------------------------------------
# treeData <- inputs$treeData
# deficit <- mastClimate( file = 'def.csv', plots = treeData$plot,
# years = treeData$year - 1, months = 6:8,
# FUN = 'sum', vname='def' )
# treeData <- cbind( treeData, deficit$x )
# summary( deficit )
## ----covars, eval=F-----------------------------------------------------------
# # include min winter temperature
# d <- "https://github.com/jimclarkatduke/mast/blob/master/tmin.csv?raw=True"
# download.file( d, destfile="tmin.csv" )
#
# # minimum winter temperature December through March of previous winter
#
# t1 <- mastClimate( file = 'tmin.csv', plots = treeData$plot,
# years = treeData$year - 1, months = 12, FUN = 'min',
# vname = 'tmin' )
# t2 <- mastClimate( file = 'tmin.csv', plots = treeData$plot,
# years = treeData$year, months = 1:3, FUN = 'min',
# vname = 'tmin' )
# tmin <- apply( cbind( t1$x[, 1], t2$x[, 1] ), 1, min )
# treeData$tminDecJanFebMar <- tmin
# inputs$treeData <- treeData
## ---- runClim, eval=F---------------------------------------------------------
# formulaRep <- as.formula( ~ diam )
# formulaFec <- as.formula( ~ diam + defJunJulAugAnom + tminDecJanFebMar )
# inputs$yearEffect <- yearEffect
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2500, burnin = 1000 )
## ----ranEff, eval=F-----------------------------------------------------------
# formulaFec <- as.formula( ~ diam ) # fecundity model
# formulaRep <- as.formula( ~ diam ) # maturation model
# inputs$randomEffect <- list( randGroups = 'tree', formulaRan = as.formula( ~ 1 ) )
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2000, burnin = 1000 )
## ----ranEff2, eval=F----------------------------------------------------------
# output <- mastif( inputs = output, ng = 4000, burnin = 1000 )
# mastPlot( output )
## ----fitSum2, eval=F----------------------------------------------------------
# output$fit
## ----regtab, eval=F-----------------------------------------------------------
# with( treeData, colSums( table( plot, region ) ) )
## ----yrpl, eval=F-------------------------------------------------------------
# yearEffect <- list( groups = c('species', 'region' ) )
## ----fit3, eval=F-------------------------------------------------------------
# inputs$treeData <- treeData
# inputs$randomEffect <- randomEffect
# inputs$yearEffect <- yearEffect
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2500, burnin = 500 )
## ----moreYR, eval=F-----------------------------------------------------------
# predList <- list( mapMeters = 10, plots = 'DUKE_BW', years = 1998:2014 )
# output$predList <- predList
# output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
# mastPlot( output )
## ----map2, eval=F-------------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/pinusExample.rdata?raw=True"
# repmis::source_data( d )
#
# mapList <- list( treeData = treeData, seedData = seedData,
# specNames = specNames, seedNames = seedNames,
# xytree = xytree, xytrap = xytrap, mapPlot = 'DUKE_EW',
# mapYears = c( 2007:2010 ), treeSymbol = treeData$diam,
# treeScale = .6, trapScale=1.4,
# plotScale = 1.2, LEGEND=T )
# mastMap( mapList )
## ----fit0, eval=F-------------------------------------------------------------
# formulaFec <- as.formula( ~ diam ) # fecundity model
# formulaRep <- as.formula( ~ diam ) # maturation model
#
# yearEffect <- list( groups = 'species', p = 4 ) # AR( 4 )
# randomEffect <- list( randGroups = 'tree',
# formulaRan = as.formula( ~ 1 ) )
#
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData,
# yearEffect = yearEffect,
# randomEffect = randomEffect,
# xytree = xytree, xytrap = xytrap, priorDist = 20,
# priorVDist = 5, minDist = 15, maxDist = 30,
# minDiam = 12, maxDiam = 40,
# maxF = 1e+6, seedTraits = seedTraits )
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 500, burnin = 100 )
## ----moreAR, eval=F-----------------------------------------------------------
# output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
# mastPlot( output, plotPars = plotPars )
## ----moreYR2, eval=F----------------------------------------------------------
# plots <- c( 'DUKE_EW', 'CWT_118' )
# years <- 1980:2025
# output$predList <- list( mapMeters = 10, plots = plots, years = years )
# output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
## ----yrPlot, eval=F-----------------------------------------------------------
# mastPlot( output, )
## ----onemap, eval=F-----------------------------------------------------------
# mapList <- output
# mapList$mapPlot <- 'DUKE_EW'
# mapList$mapYears <- c( 2011:2012 )
# mapList$PREDICT <- T
# mapList$treeScale <- .5
# mapList$trapScale <- .8
# mapList$LEGEND <- T
# mapList$scaleValue <- 50
# mapList$plotScale <- 2
# mapList$COLORSCALE <- T
# mapList$mfrow <- c( 2, 1 )
#
# mastMap( mapList )
## ----onemap1, eval=F----------------------------------------------------------
# mapList$mapPlot <- 'CWT_118'
# mapList$mapYears <- 2015
# mapList$PREDICT <- T
# mapList$treeScale <- 1.5
# mapList$trapScale <- .8
# mapList$LEGEND <- T
# mapList$scaleValue <- 50
# mapList$plotScale <- 2
# mapList$COLORSCALE <- T
# mapList$mfrow <- c( 1, 1 )
# mastMap( mapList )
## ----outpars0, eval=F---------------------------------------------------------
# summary( output )
## ---- eval=F, echo = F--------------------------------------------------------
#
# # abies output west
# load('prediction.rdata')
#
# fecPred <- prediction$fecPred
# fecPred <- fecPred[ fecPred$matrEst > .5, ]
#
# ptab <- table( fecPred$plot, fecPred$species )
# ww <- which( ptab > 1000, arr.ind = T )
# keep <- paste( rownames(ptab)[ww[,1]], colnames(ptab)[ww[,2]] )
#
# fecPred$plotSpec <- paste(fecPred$plot, fecPred$species)
# fecPred <- fecPred[ fecPred$plotSpec %in% keep, ]
# specs <- sort( unique( fecPred$species ) )
#
# yseq <- seq( 0, 10, length = 100 )
# intVal <- matrix( 0, length(specs), 100 )
# rownames( intVal ) <- specs
#
# par( mfrow = c(1, 3), bty = 'n', mar = c(2,4,1,1), omi = c(.5,.1,.1,.1) )
# plot( NA, xlim = c(2, 10), ylim = c(.01, 1), xlab = 'Year', ylab = 'Density/yr',
# log = 'xy')
#
# for( i in 1:length(keep) ){
#
# wi <- which( fecPred$plotSpec == keep[i] )
# ci <- fecPred$species[wi[1]]
# col <- match( ci, specs )
#
# tmp <- mastVolatility( treeID = fecPred$treeID[wi], year = fecPred$year[wi],
# fec = fecPred$fecEstMu[wi] )
# if( is.null(tmp) )next
#
# intVal[ col, ] <- intVal[ col, ] + dnorm( yseq, tmp$stats['Period', 1], tmp$stats['Period', 2] )
#
# dens <- tmp$statsDensity
# lines( dens[ 'Period', ], dens[ 'Mean', ], lwd = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] - dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] + dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# }
# title( 'a) Plot-species groups' )
#
# plot( NA, xlim = c(2, 6), ylim = c(0, .12), xlab = '', ylab = 'Density' )
# for( i in 1:length(specs) ){
# polygon( yseq, intVal[i,]/sum(intVal[i,]), border = i, col = getColor(i, .4) )
# }
# legend( 'topright', specs, text.col = c(1:length(specs)), bty = 'n', cex = .8 )
# title( 'b) Period estimates' )
#
# # using year effects
#
# load('parameters.rdata')
#
# # if there are year effects in the model:
#
# betaYrRand <- parameters$betaYrRand
# betaYrSE <- parameters$betaYrRandSE
#
# spec <- strsplit( rownames( betaYrRand ), '_' )
# spec <- sapply( spec, function(x) x[2] )
# specs <- sort( unique(spec) )
#
# mastMatrix <- matrix( 0, nrow(betaYrRand), 5 )
# rownames(mastMatrix) <- rownames(betaYrRand)
# colnames(mastMatrix) <- c( 'nyr', 'Variance', 'Volatility', 'Period Est', 'Period SD' )
#
#
# #par( mfrow = c(1, 2), bty = 'n' )
# plot( NA, xlim = c(2, 10), ylim = c(.002, .4), xlab = '', ylab = 'Density/yr',
# log = 'xy')
#
# for(i in 1:nrow(betaYrRand)){
#
# wc <- which( betaYrRand[i,] != 0 )
# if( length(wc) < 6 )next
#
# s <- spectralDensity( betaYrRand[i,wc] )
# if( !is.matrix( s$spect ) )next
#
# mastMatrix[i, ] <- c( length(wc), s$totVar, s$volatility, s$periodMu, s$periodSd )
#
# period <- 1/s$spec[, 'frequency' ]
# dens <- s$spec[, 'spectralDensity' ]/length(wc) # series vary in length
#
# col <- match( spec[i], specs )
#
# lines( period, dens, lwd = 2, col = getColor( col, .4) )
# }
# title( 'c) Year effects' )
# mtext( 'Period (yr)', 1, line = 1, outer = T )
#
# keepRows <- which( is.finite(mastMatrix[,'Variance']) & mastMatrix[,'Variance'] != 0 )
# keepCols <- which( colSums( frequency, na.rm=T ) > 0 )
# mastMatrix <- mastMatrix[ keepRows, ]
#
# save( fecPred, betaYrRand, file = 'outputAbies.rdata' )
## ---- eval = F----------------------------------------------------------------
# getColor <- function( col, trans ){ # transparent colors
# tmp <- col2rgb( col )
# rgb( tmp[ 1, ], tmp[ 2, ], tmp[ 3, ], maxColorValue = 255,
# alpha = 255*trans, names = paste( col, trans, sep = '_' ) )
# }
#
# # from output$prediction$fecPred and output$parameters$betaYrRand
# d <- "https://github.com/jimclarkatduke/mast/blob/master/outputAbies.rdata?raw=True"
# repmis::source_data( d )
#
# specs <- sort( unique( fecPred$species ) ) # accumulate period estimates
# yseq <- seq( 0, 10, length = 100 )
# intVal <- matrix( 0, length(specs), 100 )
# rownames( intVal ) <- specs
# plotSpecs <- sort( unique( fecPred$plotSpec ) ) # label trees in a plot-species group
#
# par( mfrow = c(1, 3), bty = 'n', mar = c(2,4,1,1), omi = c(.5,.1,.1,.1) )
#
# plot( NA, xlim = c(2, 20), ylim = c(.01, 1), xlab = '', ylab = 'Density/nyr', log = 'xy')
#
# for( i in 1:length(plotSpecs) ){
#
# wi <- which( fecPred$plotSpec == plotSpecs[i] ) # tree-years in group
# ci <- fecPred$species[wi[1]]
# col <- match( ci, specs ) # color by species
#
# tmp <- mastVolatility( treeID = fecPred$treeID[wi], year = fecPred$year[wi],
# fec = fecPred$fecEstMu[wi] )
# if( is.null(tmp) )next
#
# intVal[ col, ] <- intVal[ col, ] + dnorm( yseq, tmp$stats['Period', 1], tmp$stats['Period', 2] )
#
# # density +/- 1 SD
# dens <- tmp$statsDensity
# lines( dens[ 'Period', ], dens[ 'Mean', ], lwd = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] - dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] + dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# }
# title( 'a) Plot-species groups' )
# legend( 'topright', specs, text.col = c(1:length(specs)), bty = 'n', cex = .8 )
#
# plot( NA, xlim = c(2, 8), ylim = c(0, .12), xlab = '', ylab = 'Density' ) # density of mean intervals
# for( i in 1:length(specs) ){
# polygon( yseq, intVal[i,]/sum(intVal[i,]), border = i, col = getColor(i, .4) )
# }
# title( 'b) Period estimates' )
## ---- eval = F----------------------------------------------------------------
# spec <- strsplit( rownames( betaYrRand ), '_' ) # extract species from ecoregion_species rownames
# spec <- sapply( spec, function(x) x[2] )
# specs <- sort( unique(spec) )
#
# mastMatrix <- matrix( 0, nrow(betaYrRand), 5 ) # store stats by group
# rownames(mastMatrix) <- rownames(betaYrRand)
# colnames(mastMatrix) <- c( 'nyr', 'Variance', 'Volatility', 'Period Est', 'Period SD' )
#
# plot( NA, xlim = c(2, 10), ylim = c(.002, .4), xlab = '', ylab = 'Density/yr', log = 'xy')
#
# for(i in 1:nrow(betaYrRand)){
#
# wc <- which( betaYrRand[i,] != 0 )
# if( length(wc) < 6 )next
#
# s <- mastSpectralDensity( betaYrRand[i,wc] )
# if( !is.matrix( s$spect ) )next
#
# mastMatrix[i, ] <- c( length(wc), s$totVar, s$volatility, s$periodMu, s$periodSd )
#
# period <- 1/s$spec[, 'frequency' ]
# dens <- s$spec[, 'spectralDensity' ]/length(wc) # series vary in length
# col <- match( spec[i], specs )
#
# lines( period, dens, lwd = 2, col = getColor( col, .4) )
# }
# title( 'c) Year effects' )
# mtext( 'Period (yr)', 1, line = 1, outer = T )
#
# keepRows <- which( is.finite(mastMatrix[,'Variance']) & mastMatrix[,'Variance'] != 0 )
# keepCols <- which( colSums( frequency, na.rm=T ) > 0 )
# mastMatrix <- mastMatrix[ keepRows, ]
Try the mastif package in your browser
Any scripts or data that you put into this service are public.
mastif documentation built on Feb. 16, 2023, 5:30 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----outform, echo=F----------------------------------------------------------
insertPlot <- function( file, caption ){
# outputFormat = knitr::opts_knit$get( "rmarkdown.pandoc.to" )
# if( outputFormat == 'latex' )
# paste( "![ ", caption, " ]( ", file, " )", sep="" )
}
bigskip <- function( ){
# outputFormat = knitr::opts_knit$get( "rmarkdown.pandoc.to" )
# if( outputFormat == 'latex' )
# "\\bigskip"
# else
"<br>"
}
## ----getFiles, eval = F, echo=F-----------------------------------------------
# Rcpp::sourceCpp( '../RcppFunctions/cppFns.cpp' )
# source( '../RFunctions/mastifFunctions.r' )
# library( RANN )
# library( robustbase )
## ----simSetup0, eval = F------------------------------------------------------
# seedNames <- specNames <- 'acerRubr'
# sim <- list( nplot=5, nyr=10, ntree=30, ntrap=40, specNames = specNames, seedNames = seedNames )
## ----sim0, eval = F-----------------------------------------------------------
# inputs <- mastSim( sim ) # simulate dispersal data
# seedData <- inputs$seedData # year, plot, trap, seed counts
# treeData <- inputs$treeData # year, plot, tree data
# xytree <- inputs$xytree # tree locations
# xytrap <- inputs$xytrap # trap locations
# formulaFec <- inputs$formulaFec # fecundity model
# formulaRep <- inputs$formulaRep # maturation model
# trueValues <- inputs$trueValues # true states and parameter values
## ----formulaFec, eval = F-----------------------------------------------------
# formulaFec
## ----treeData0, eval = F------------------------------------------------------
# head( treeData )
## ----xytree, eval = F---------------------------------------------------------
# head( xytree, 5 )
## ----treeData1, eval = F------------------------------------------------------
# head( seedData )
## ----map1a, eval = F----------------------------------------------------------
# dataTab <- table( treeData$plot, treeData$year )
#
# w <- which( dataTab > 0, arr.ind=T ) # a plot-year with observations
# w <- w[sample( nrow( w ), 1 ), ]
#
# mapYears <- as.numeric( colnames( dataTab )[w[2]] )
# mapPlot <- rownames( dataTab )[w[1]]
# inputs$mapPlot <- mapPlot
# inputs$mapYears <- mapYears
# inputs$treeSymbol <- treeData$diam
# inputs$SCALEBAR <- T
#
# mastMap( inputs )
## ----map3, eval=F-------------------------------------------------------------
# inputs$treeSymbol <- trueValues$fec
# inputs$treeScale <- 2
# inputs$trapScale <- 1
# mastMap( inputs )
## ----map4, eval=F-------------------------------------------------------------
# inputs$treeSymbol <- trueValues$repr
# inputs$treeScale <- .5
# mastMap( inputs )
## ----hist0, eval = F----------------------------------------------------------
# par( mfrow=c( 2, 1 ), bty='n', mar=c( 4, 4, 1, 1 ) )
# seedData <- inputs$seedData
# seedNames <- inputs$seedNames
#
# hist( as.matrix( seedData[, seedNames] ) , nclass=100,
# xlab = 'seed count', ylab='per trap', main='' )
# hist( trueValues$fec, nclass=100, xlab = 'seeds produced', ylab = 'per tree', main = '' )
## ----mast2, eval = F----------------------------------------------------------
# output <- mastif( inputs = inputs, ng = 4000, burnin = 500 )
## ----tabPars0, eval = F-------------------------------------------------------
# summary( output )
## ----pars, eval = F-----------------------------------------------------------
# plotPars <- list( trueValues = trueValues )
# mastPlot( output, plotPars )
## ---- eval=F------------------------------------------------------------------
# output <- mastif( inputs = output, ng = 5000, burnin = 3000 )
## ----restart, eval=F----------------------------------------------------------
# inputs$predList <- list( mapMeters = 3, plots = mapPlot, years = mapYears )
# output <- mastif( inputs = output, ng = 2000, burnin = 1000 )
## ----mapout, eval = F---------------------------------------------------------
# output$mapPlot <- mapPlot
# output$mapYears <- mapYears
# output$treeScale <- 1.5
# output$trapScale <- .8
# output$PREDICT <- T
# output$scaleValue <- 20
# output$plotScale <- 1
# output$COLORSCALE <- T
# output$LEGEND <- T
#
# mastMap( output )
## ----to rmd, eval=F-----------------------------------------------------------
# plotPars$RMD <- 'pdf'
# mastPlot( output, plotPars )
## ----simSetup, eval = F-------------------------------------------------------
# specNames <- c( 'pinuTaeda', 'pinuEchi', 'pinuVirg' )
# seedNames <- c( 'pinuTaeda', 'pinuEchi', 'pinuVirg', 'pinuUNKN' )
# sim <- list( nyr=4, ntree=25, nplot=10, ntrap=50, specNames = specNames,
# seedNames = seedNames )
## ----sim, eval = F------------------------------------------------------------
# inputs <- mastSim( sim ) # simulate dispersal data
# R <- inputs$trueValues$R # species to seedNames probability matrix
# round( R, 2 )
## ----mast3, eval = F----------------------------------------------------------
# output <- mastif( inputs = inputs, ng = 2000, burnin = 1000 )
## ----tabPars, eval = F--------------------------------------------------------
# summary( output )
## ----pars0, eval = F----------------------------------------------------------
# plotPars <- list( trueValues = inputs$trueValues )
# mastPlot( output, plotPars )
## ----again, eval=F------------------------------------------------------------
# tab <- with( inputs$seedData, table( plot, year ) )
# years <- as.numeric( colnames( tab )[tab[1, ] > 0] ) # years for 1st plot
# output$predList <- list( plots = 'p1', years = years )
# output <- mastif( inputs = output, ng = 3000, burnin = 1500 )
# mastPlot( output )
## ---- eval=F------------------------------------------------------------------
# specNames <- c( 'pinuTaed', 'pinuEchi' )
#
# #seeds never differentiated:
# seedNames <- c( 'pinuUNKN' )
#
# #one species sometimes differentiated:
# seedNames <- c( 'pinuTaed', 'pinuUNKN' )
#
# #both species sometimes differentiated:
# seedNames <- c( 'pinuTaed', 'pinuEchi', 'pinuUNKN' )
## ----map21, eval=F------------------------------------------------------------
# library( repmis )
# d <- "https://github.com/jimclarkatduke/mast/blob/master/liriodendronExample.rData?raw=True"
# repmis::source_data( d )
# mapList <- list( treeData = treeData, seedData = seedData,
# specNames = specNames, seedNames = seedNames,
# xytree = xytree, xytrap = xytrap, mapPlot = 'DUKE_BW',
# mapYears = 2011:2014, treeSymbol = treeData$diam,
# treeScale = .7, trapScale = 1.5, plotScale = 2,
# SCALEBAR=T, scaleValue=50 )
# mastMap( mapList )
## ----litu1, eval=F------------------------------------------------------------
# head( treeData, 3 )
## ----litu2, eval=F------------------------------------------------------------
# head( seedData, 3 )
## ----fit, eval=F--------------------------------------------------------------
# formulaFec <- as.formula( ~ diam ) # fecundity model
# formulaRep <- as.formula( ~ diam ) # maturation model
#
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData, xytree = xytree,
# xytrap = xytrap )
# output <- mastif( inputs, formulaFec, formulaRep, ng = 3000, burnin = 1000 )
## ----more, eval=F-------------------------------------------------------------
# output$predList <- list( mapMeters = 10, plots = 'DUKE_EW', years = 2010:2015 )
# output <- mastif( inputs = output, ng = 4000, burnin = 1000 )
## ----plotmydata1, eval=F------------------------------------------------------
# mastPlot( output )
## ----outpars, eval=F----------------------------------------------------------
# summary( output )
## ----fitSum, eval=F-----------------------------------------------------------
# output$fit
## ----nopred, eval=F-----------------------------------------------------------
# #group plots in regions for year effects
# region <- rep( 'sApps', nrow( treeData ) )
# region[ as.character( treeData$plot ) == 'DUKE_EW' ] <- 'piedmont'
#
# treeData$region <- region
#
# formulaFec <- as.formula( ~ diam )
# formulaRep <- as.formula( ~ diam )
# yearEffect <- list( groups = 'region' )
# randomEffect <- list( randGroups = 'treeID', formulaRan = as.formula( ~ 1 ) )
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData,
# xytree = xytree, xytrap = xytrap,
# priorDist = 28, priorVDist = 15, maxDist = 50, minDist = 15,
# minDiam = 25, maxF = 1e+6,
# randomEffect = randomEffect, yearEffect = yearEffect )
# output <- mastif( inputs, formulaFec, formulaRep, ng = 2000, burnin = 1000 )
# mastPlot( output )
## ---- yr, eval=F--------------------------------------------------------------
# yearEffect <- list( groups = c( 'species', 'region' ) ) # year effects
## ---- eval=F------------------------------------------------------------------
# inputs$priorList <- list( minDiam = 15, maxDiam = 60 )
## ---- eval=F------------------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/priorParametersPinus.csv?raw=True"
# download.file( d, destfile="priorParametersPinus.csv" )
#
# specNames <- c( "pinuEchi", "pinuRigi", "pinuStro", "pinuTaed", "pinuVirg" )
# seedNames <- c( "pinuEchi", "pinuRigi", "pinuStro", "pinuTaed", "pinuVirg", "pinuUNKN" )
# priorTable <- mastPriors( "priorParametersPinus.csv", specNames,
# code = 'code4', genus = 'pinus' )
## ---- eval=F, echo=F----------------------------------------------------------
#
# # THIS BLOCK IS OMITTED
#
# load( '../combinedSites/pinus.Rdata' )
#
# pl <- c( 'CWT_118', 'CWT_218', 'DUKE_BW', 'DUKE_EW', 'MARS_F', 'MARS_P' )
# treeData <- treeData[treeData$plot %in% pl, ]
# seedData <- seedData[seedData$plot %in% pl, ]
# xytree <- xytree[xytree$plot %in% pl, ]
# xytrap <- xytrap[xytrap$plot %in% pl, ]
#
# count <- as.matrix( seedData[, -c( 1:5 )] )
# count <- count[, colSums( count, na.rm=T ) > 0]
# count <- count[, -1]
# seedData <- cbind( seedData[, 1:5], count )
#
# treeData <- treeData[, 1:6]
#
# seedNames <- colnames( count )
# specNames <- sort( unique( treeData$species ) )
#
# save( treeData, seedData, xytree, xytrap,
# specNames, seedNames, file='pinusExample.rData' )
#
# #load( 'pinusExample.rData' )
## ----priorBeta, eval=F--------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/pinusExample.rdata?raw=True"
# repmis::source_data( d )
## ---- eval=F------------------------------------------------------------------
# formulaRep <- as.formula( ~ diam )
# formulaFec <- as.formula( ~ diam + I( diam^2 ) )
# betaPrior <- list( pos = 'diam', neg = 'I( diam^2 )' )
#
# inputs <- list( treeData = treeData, seedData = seedData, xytree = xytree,
# xytrap = xytrap, specNames = specNames, seedNames = seedNames,
# betaPrior = betaPrior, priorTable = priorTable,
# seedTraits = seedTraits )
# output <- mastif( inputs = inputs, formulaFec, formulaRep,
# ng = 500, burnin = 200 )
#
# # restart
# output <- mastif( output, formulaFec, formulaRep, ng = 5000, burnin = 1000 )
# mastPlot( output )
## ----fill, eval=F-------------------------------------------------------------
# # randomly remove years for this example:
# years <- sort( unique( treeData$year ) )
# sy <- sample( years, 5 )
# treeData <- treeData[treeData$year %in% sy, ]
# treeData[1:10, ]
## ---- removed, eval=F---------------------------------------------------------
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData,
# xytree = xytree, xytrap = xytrap, priorTable = priorTable,
# seedTraits = seedTraits )
# inputs <- mastFillCensus( inputs, beforeFirst=10, afterLast=10 )
# inputs$treeData[1:10, ]
## ---- treeYr table, eval=F----------------------------------------------------
# # original data
# table( treeData$year )
#
# # filled census
# table( inputs$treeData$year )
# table( seedData$year )
## ---- treeOnly, eval=F--------------------------------------------------------
# p <- 3
# inputs <- mastFillCensus( inputs, p = p )
# treeData <- inputs$treeData
## ----arr, eval=F--------------------------------------------------------------
# region <- c( 'CWT', 'DUKE', 'HARV' )
# treeData$region <- 'SCBI'
# for( j in 1:length( region ) ){
# wj <- which( startsWith( treeData$plot, region[j] ) )
# treeData$region[wj] <- region[j]
# }
# inputs$treeData <- treeData
# yearEffect <- list( groups = c( 'species', 'region' ), p = p )
## ----def, eval=F--------------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/def.csv?raw=True"
# download.file( d, destfile="def.csv" )
## ----types, eval=F------------------------------------------------------------
# treeData <- inputs$treeData
# deficit <- mastClimate( file = 'def.csv', plots = treeData$plot,
# years = treeData$year - 1, months = 6:8,
# FUN = 'sum', vname='def' )
# treeData <- cbind( treeData, deficit$x )
# summary( deficit )
## ----covars, eval=F-----------------------------------------------------------
# # include min winter temperature
# d <- "https://github.com/jimclarkatduke/mast/blob/master/tmin.csv?raw=True"
# download.file( d, destfile="tmin.csv" )
#
# # minimum winter temperature December through March of previous winter
#
# t1 <- mastClimate( file = 'tmin.csv', plots = treeData$plot,
# years = treeData$year - 1, months = 12, FUN = 'min',
# vname = 'tmin' )
# t2 <- mastClimate( file = 'tmin.csv', plots = treeData$plot,
# years = treeData$year, months = 1:3, FUN = 'min',
# vname = 'tmin' )
# tmin <- apply( cbind( t1$x[, 1], t2$x[, 1] ), 1, min )
# treeData$tminDecJanFebMar <- tmin
# inputs$treeData <- treeData
## ---- runClim, eval=F---------------------------------------------------------
# formulaRep <- as.formula( ~ diam )
# formulaFec <- as.formula( ~ diam + defJunJulAugAnom + tminDecJanFebMar )
# inputs$yearEffect <- yearEffect
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2500, burnin = 1000 )
## ----ranEff, eval=F-----------------------------------------------------------
# formulaFec <- as.formula( ~ diam ) # fecundity model
# formulaRep <- as.formula( ~ diam ) # maturation model
# inputs$randomEffect <- list( randGroups = 'tree', formulaRan = as.formula( ~ 1 ) )
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2000, burnin = 1000 )
## ----ranEff2, eval=F----------------------------------------------------------
# output <- mastif( inputs = output, ng = 4000, burnin = 1000 )
# mastPlot( output )
## ----fitSum2, eval=F----------------------------------------------------------
# output$fit
## ----regtab, eval=F-----------------------------------------------------------
# with( treeData, colSums( table( plot, region ) ) )
## ----yrpl, eval=F-------------------------------------------------------------
# yearEffect <- list( groups = c('species', 'region' ) )
## ----fit3, eval=F-------------------------------------------------------------
# inputs$treeData <- treeData
# inputs$randomEffect <- randomEffect
# inputs$yearEffect <- yearEffect
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 2500, burnin = 500 )
## ----moreYR, eval=F-----------------------------------------------------------
# predList <- list( mapMeters = 10, plots = 'DUKE_BW', years = 1998:2014 )
# output$predList <- predList
# output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
# mastPlot( output )
## ----map2, eval=F-------------------------------------------------------------
# d <- "https://github.com/jimclarkatduke/mast/blob/master/pinusExample.rdata?raw=True"
# repmis::source_data( d )
#
# mapList <- list( treeData = treeData, seedData = seedData,
# specNames = specNames, seedNames = seedNames,
# xytree = xytree, xytrap = xytrap, mapPlot = 'DUKE_EW',
# mapYears = c( 2007:2010 ), treeSymbol = treeData$diam,
# treeScale = .6, trapScale=1.4,
# plotScale = 1.2, LEGEND=T )
# mastMap( mapList )
## ----fit0, eval=F-------------------------------------------------------------
# formulaFec <- as.formula( ~ diam ) # fecundity model
# formulaRep <- as.formula( ~ diam ) # maturation model
#
# yearEffect <- list( groups = 'species', p = 4 ) # AR( 4 )
# randomEffect <- list( randGroups = 'tree',
# formulaRan = as.formula( ~ 1 ) )
#
# inputs <- list( specNames = specNames, seedNames = seedNames,
# treeData = treeData, seedData = seedData,
# yearEffect = yearEffect,
# randomEffect = randomEffect,
# xytree = xytree, xytrap = xytrap, priorDist = 20,
# priorVDist = 5, minDist = 15, maxDist = 30,
# minDiam = 12, maxDiam = 40,
# maxF = 1e+6, seedTraits = seedTraits )
# output <- mastif( inputs = inputs, formulaFec, formulaRep, ng = 500, burnin = 100 )
## ----moreAR, eval=F-----------------------------------------------------------
# output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
# mastPlot( output, plotPars = plotPars )
## ----moreYR2, eval=F----------------------------------------------------------
# plots <- c( 'DUKE_EW', 'CWT_118' )
# years <- 1980:2025
# output$predList <- list( mapMeters = 10, plots = plots, years = years )
# output <- mastif( inputs = output, ng = 3000, burnin = 1000 )
## ----yrPlot, eval=F-----------------------------------------------------------
# mastPlot( output, )
## ----onemap, eval=F-----------------------------------------------------------
# mapList <- output
# mapList$mapPlot <- 'DUKE_EW'
# mapList$mapYears <- c( 2011:2012 )
# mapList$PREDICT <- T
# mapList$treeScale <- .5
# mapList$trapScale <- .8
# mapList$LEGEND <- T
# mapList$scaleValue <- 50
# mapList$plotScale <- 2
# mapList$COLORSCALE <- T
# mapList$mfrow <- c( 2, 1 )
#
# mastMap( mapList )
## ----onemap1, eval=F----------------------------------------------------------
# mapList$mapPlot <- 'CWT_118'
# mapList$mapYears <- 2015
# mapList$PREDICT <- T
# mapList$treeScale <- 1.5
# mapList$trapScale <- .8
# mapList$LEGEND <- T
# mapList$scaleValue <- 50
# mapList$plotScale <- 2
# mapList$COLORSCALE <- T
# mapList$mfrow <- c( 1, 1 )
# mastMap( mapList )
## ----outpars0, eval=F---------------------------------------------------------
# summary( output )
## ---- eval=F, echo = F--------------------------------------------------------
#
# # abies output west
# load('prediction.rdata')
#
# fecPred <- prediction$fecPred
# fecPred <- fecPred[ fecPred$matrEst > .5, ]
#
# ptab <- table( fecPred$plot, fecPred$species )
# ww <- which( ptab > 1000, arr.ind = T )
# keep <- paste( rownames(ptab)[ww[,1]], colnames(ptab)[ww[,2]] )
#
# fecPred$plotSpec <- paste(fecPred$plot, fecPred$species)
# fecPred <- fecPred[ fecPred$plotSpec %in% keep, ]
# specs <- sort( unique( fecPred$species ) )
#
# yseq <- seq( 0, 10, length = 100 )
# intVal <- matrix( 0, length(specs), 100 )
# rownames( intVal ) <- specs
#
# par( mfrow = c(1, 3), bty = 'n', mar = c(2,4,1,1), omi = c(.5,.1,.1,.1) )
# plot( NA, xlim = c(2, 10), ylim = c(.01, 1), xlab = 'Year', ylab = 'Density/yr',
# log = 'xy')
#
# for( i in 1:length(keep) ){
#
# wi <- which( fecPred$plotSpec == keep[i] )
# ci <- fecPred$species[wi[1]]
# col <- match( ci, specs )
#
# tmp <- mastVolatility( treeID = fecPred$treeID[wi], year = fecPred$year[wi],
# fec = fecPred$fecEstMu[wi] )
# if( is.null(tmp) )next
#
# intVal[ col, ] <- intVal[ col, ] + dnorm( yseq, tmp$stats['Period', 1], tmp$stats['Period', 2] )
#
# dens <- tmp$statsDensity
# lines( dens[ 'Period', ], dens[ 'Mean', ], lwd = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] - dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] + dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# }
# title( 'a) Plot-species groups' )
#
# plot( NA, xlim = c(2, 6), ylim = c(0, .12), xlab = '', ylab = 'Density' )
# for( i in 1:length(specs) ){
# polygon( yseq, intVal[i,]/sum(intVal[i,]), border = i, col = getColor(i, .4) )
# }
# legend( 'topright', specs, text.col = c(1:length(specs)), bty = 'n', cex = .8 )
# title( 'b) Period estimates' )
#
# # using year effects
#
# load('parameters.rdata')
#
# # if there are year effects in the model:
#
# betaYrRand <- parameters$betaYrRand
# betaYrSE <- parameters$betaYrRandSE
#
# spec <- strsplit( rownames( betaYrRand ), '_' )
# spec <- sapply( spec, function(x) x[2] )
# specs <- sort( unique(spec) )
#
# mastMatrix <- matrix( 0, nrow(betaYrRand), 5 )
# rownames(mastMatrix) <- rownames(betaYrRand)
# colnames(mastMatrix) <- c( 'nyr', 'Variance', 'Volatility', 'Period Est', 'Period SD' )
#
#
# #par( mfrow = c(1, 2), bty = 'n' )
# plot( NA, xlim = c(2, 10), ylim = c(.002, .4), xlab = '', ylab = 'Density/yr',
# log = 'xy')
#
# for(i in 1:nrow(betaYrRand)){
#
# wc <- which( betaYrRand[i,] != 0 )
# if( length(wc) < 6 )next
#
# s <- spectralDensity( betaYrRand[i,wc] )
# if( !is.matrix( s$spect ) )next
#
# mastMatrix[i, ] <- c( length(wc), s$totVar, s$volatility, s$periodMu, s$periodSd )
#
# period <- 1/s$spec[, 'frequency' ]
# dens <- s$spec[, 'spectralDensity' ]/length(wc) # series vary in length
#
# col <- match( spec[i], specs )
#
# lines( period, dens, lwd = 2, col = getColor( col, .4) )
# }
# title( 'c) Year effects' )
# mtext( 'Period (yr)', 1, line = 1, outer = T )
#
# keepRows <- which( is.finite(mastMatrix[,'Variance']) & mastMatrix[,'Variance'] != 0 )
# keepCols <- which( colSums( frequency, na.rm=T ) > 0 )
# mastMatrix <- mastMatrix[ keepRows, ]
#
# save( fecPred, betaYrRand, file = 'outputAbies.rdata' )
## ---- eval = F----------------------------------------------------------------
# getColor <- function( col, trans ){ # transparent colors
# tmp <- col2rgb( col )
# rgb( tmp[ 1, ], tmp[ 2, ], tmp[ 3, ], maxColorValue = 255,
# alpha = 255*trans, names = paste( col, trans, sep = '_' ) )
# }
#
# # from output$prediction$fecPred and output$parameters$betaYrRand
# d <- "https://github.com/jimclarkatduke/mast/blob/master/outputAbies.rdata?raw=True"
# repmis::source_data( d )
#
# specs <- sort( unique( fecPred$species ) ) # accumulate period estimates
# yseq <- seq( 0, 10, length = 100 )
# intVal <- matrix( 0, length(specs), 100 )
# rownames( intVal ) <- specs
# plotSpecs <- sort( unique( fecPred$plotSpec ) ) # label trees in a plot-species group
#
# par( mfrow = c(1, 3), bty = 'n', mar = c(2,4,1,1), omi = c(.5,.1,.1,.1) )
#
# plot( NA, xlim = c(2, 20), ylim = c(.01, 1), xlab = '', ylab = 'Density/nyr', log = 'xy')
#
# for( i in 1:length(plotSpecs) ){
#
# wi <- which( fecPred$plotSpec == plotSpecs[i] ) # tree-years in group
# ci <- fecPred$species[wi[1]]
# col <- match( ci, specs ) # color by species
#
# tmp <- mastVolatility( treeID = fecPred$treeID[wi], year = fecPred$year[wi],
# fec = fecPred$fecEstMu[wi] )
# if( is.null(tmp) )next
#
# intVal[ col, ] <- intVal[ col, ] + dnorm( yseq, tmp$stats['Period', 1], tmp$stats['Period', 2] )
#
# # density +/- 1 SD
# dens <- tmp$statsDensity
# lines( dens[ 'Period', ], dens[ 'Mean', ], lwd = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] - dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# lines( dens[ 'Period', ], dens[ 'Mean', ] + dens[ 'SD', ], lty = 2, col = getColor( col, .4) )
# }
# title( 'a) Plot-species groups' )
# legend( 'topright', specs, text.col = c(1:length(specs)), bty = 'n', cex = .8 )
#
# plot( NA, xlim = c(2, 8), ylim = c(0, .12), xlab = '', ylab = 'Density' ) # density of mean intervals
# for( i in 1:length(specs) ){
# polygon( yseq, intVal[i,]/sum(intVal[i,]), border = i, col = getColor(i, .4) )
# }
# title( 'b) Period estimates' )
## ---- eval = F----------------------------------------------------------------
# spec <- strsplit( rownames( betaYrRand ), '_' ) # extract species from ecoregion_species rownames
# spec <- sapply( spec, function(x) x[2] )
# specs <- sort( unique(spec) )
#
# mastMatrix <- matrix( 0, nrow(betaYrRand), 5 ) # store stats by group
# rownames(mastMatrix) <- rownames(betaYrRand)
# colnames(mastMatrix) <- c( 'nyr', 'Variance', 'Volatility', 'Period Est', 'Period SD' )
#
# plot( NA, xlim = c(2, 10), ylim = c(.002, .4), xlab = '', ylab = 'Density/yr', log = 'xy')
#
# for(i in 1:nrow(betaYrRand)){
#
# wc <- which( betaYrRand[i,] != 0 )
# if( length(wc) < 6 )next
#
# s <- mastSpectralDensity( betaYrRand[i,wc] )
# if( !is.matrix( s$spect ) )next
#
# mastMatrix[i, ] <- c( length(wc), s$totVar, s$volatility, s$periodMu, s$periodSd )
#
# period <- 1/s$spec[, 'frequency' ]
# dens <- s$spec[, 'spectralDensity' ]/length(wc) # series vary in length
# col <- match( spec[i], specs )
#
# lines( period, dens, lwd = 2, col = getColor( col, .4) )
# }
# title( 'c) Year effects' )
# mtext( 'Period (yr)', 1, line = 1, outer = T )
#
# keepRows <- which( is.finite(mastMatrix[,'Variance']) & mastMatrix[,'Variance'] != 0 )
# keepCols <- which( colSums( frequency, na.rm=T ) > 0 )
# mastMatrix <- mastMatrix[ keepRows, ]
Try the mastif package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.