text:
require(randomForest)
require(caret)

# Create a clean random forest ready data frame
createRFData <-
function(data, outcomeColName, numericPredictors, categoricalPredictors) {
    rfData <- data[,which(names(data) %in% numericPredictors)]
    for (i in 1:length(categoricalPredictors)) {
        rfData[, categoricalPredictors[i]] <- as.factor(data[, categoricalPredictors[i]])
    }
    rfData[, outcomeColName] <- as.factor(data[, outcomeColName])
    return(rfData)
}

# Function to find the important variables based on the mini Gini decrease
findImportantVariables <-
function(data, outcomeColName) {
    colNames <- names(data)
    outcomeCol <- which(colNames == outcomeColName)

    model <- randomForest(
        data[, outcomeCol] ~ .,
        data=data[, -outcomeCol],
        ntree=10)

    imp <- importance(initialModel)
    imp <- imp[order(-imp),]
    impDf <- data.frame(MeanDecreaseGini=imp)
    impDf$VariableName <- factor(names(imp), levels <- rev(names(imp)))

    # varImpPlot(initialModel) produces the plot but I want prettier one:
    # impPlot <- ggplot(data=impDf, aes(x=VariableName, y=MeanDecreaseGini)) +
    #    geom_bar(fill="salmon", stat="identity") + coord_flip()

    # Also create log version to show where to cut off
    # logImpPlot <- ggplot(data=impDf, aes(x=VariableName, y=log(MeanDecreaseGini))) +
    #     geom_point(fill="salmon", stat="identity") +
    #     stat_smooth(
    #         method <- "lm",
    #         data=impDf,
    #         aes(x=as.numeric(VariableName),
    #             y=log(MeanDecreaseGini))) +
    #     coord_flip()

    # Based on the mean Gini decrease analysis, use the most important predictors
    logImpDf <- data.frame(x=as.numeric(impDf$VariableName), y=log(impDf$MeanDecreaseGini))
    lmLogImp <- lm(y ~ x, data=logImpDf)
    newData <- data.frame(x=c(nrow(logImpDf):1))
    predictedLogImp <- predict(lmLogImp, newdata=newData)

    # Count until the predicted log(MeanDecreaseGini) becomes greater than actual
    numImp <- 0
    for (i in 1:length(newData$x)) {
        if (predictedLogImp[i] > logImpDf[i,]$y) break
        numImp <- numImp + 1
    }
    numImp <- length(newData$x)

    # Use top numImp predictors for this model
    predictorsUsed <- names(imp)[1:numImp]
    return (prdictorsUsed)
}

# Do K-fold cross-validation and return the confusion matrix
crossValidationConfusionMatrix <-
function(data, outcomeColName, k=10) {
    colNames <- names(data)
    outcomeCol <- which(colNames == outcomeColName)

    # Assign each observation to one of 10 folds
    # Note that this replace=TRUE is NOT replacing the observation when sampling
    id <- sample(1:k, nrow(data), replace=TRUE)
    list <- 1:k

    prediction <- data.frame()
    actual <- data.frame()

    for (i in 1:k){
        # Create training set from all training date except i-th fold
        currentTraining <- subset(data, id %in% list[-i])
        currentTesting <- subset(data, id %in% c(i))

        currentModel <- randomForest(
            currentTraining[, outcomeCol] ~ .,
            data=currentTraining[, -outcomeCol],
            ntree=10)

        message(paste("Model ", i))
        # print(currentModel$confusion)

        currentPrediction <- as.data.frame(predict(currentModel, currentTesting))
        prediction <- rbind(prediction, currentPrediction)

        currentActual <- as.data.frame(currentTesting[, outcomeCol])
        actual <- rbind(actual, currentActual)
    }

    result <- cbind(actual[, 1], prediction)
    names(result) <- c("Actual", "Predicted")
    confusionMatrix <- table(result)
    return(confusionMatrix)
}

max_stars_repo_name: daigotanaka/r-utils

max_issues_repo_name: daigotanaka/r-utils

max_forks_repo_name: daigotanaka/r-utils


text:
# Looker API 3.0 Reference
#
# ### Authorization  The Looker API uses Looker **API3** credentials for authorization and access control. Looker admins can create API3 credentials on Looker's **Admin/Users** page. Pass API3 credentials to the **/login** endpoint to obtain a temporary access_token. Include that access_token in the Authorization header of Looker API requests. For details, see [Looker API Authorization](https://looker.com/docs/r/api/authorization)  ### Client SDKs  The Looker API is a RESTful system that should be usable by any programming language capable of making HTTPS requests. Client SDKs for a variety of programming languages can be generated from the Looker API's Swagger JSON metadata to streamline use of the Looker API in your applications. A client SDK for Ruby is available as an example. For more information, see [Looker API Client SDKs](https://looker.com/docs/r/api/client_sdks)  ### Try It Out!  The 'api-docs' page served by the Looker instance includes 'Try It Out!' buttons for each API method. After logging in with API3 credentials, you can use the \"Try It Out!\" buttons to call the API directly from the documentation page to interactively explore API features and responses.  ### Versioning  Future releases of Looker will expand this API release-by-release to securely expose more and more of the core power of Looker to API client applications. API endpoints marked as \"beta\" may receive breaking changes without warning. Stable (non-beta) API endpoints should not receive breaking changes in future releases. For more information, see [Looker API Versioning](https://looker.com/docs/r/api/versioning)
#
# OpenAPI spec version: 3.0.0
#
# Generated by: https://github.com/swagger-api/swagger-codegen.git

#' ContentMeta Class
#'
#' @field id
#' @field name
#' @field parent_id
#' @field space_id
#' @field dashboard_id
#' @field look_id
#' @field content_type
#' @field inherits
#' @field inheriting_id
#' @field can
#'
#' @importFrom R6 R6Class
#' @importFrom jsonlite fromJSON toJSON
#' @export
ContentMeta <- R6::R6Class(
  'ContentMeta',
  public = list(
    `id` = NULL,
    `name` = NULL,
    `parent_id` = NULL,
    `space_id` = NULL,
    `dashboard_id` = NULL,
    `look_id` = NULL,
    `content_type` = NULL,
    `inherits` = NULL,
    `inheriting_id` = NULL,
    `can` = NULL,
    initialize = function(`id`, `name`, `parent_id`, `space_id`, `dashboard_id`, `look_id`, `content_type`, `inherits`, `inheriting_id`, `can`){
      if (!missing(`id`)) {
        stopifnot(is.numeric(`id`), length(`id`) == 1)
        self$`id` <- `id`
      }
      if (!missing(`name`)) {
        stopifnot(is.character(`name`), length(`name`) == 1)
        self$`name` <- `name`
      }
      if (!missing(`parent_id`)) {
        stopifnot(is.numeric(`parent_id`), length(`parent_id`) == 1)
        self$`parent_id` <- `parent_id`
      }
      if (!missing(`space_id`)) {
        stopifnot(is.character(`space_id`), length(`space_id`) == 1)
        self$`space_id` <- `space_id`
      }
      if (!missing(`dashboard_id`)) {
        stopifnot(is.character(`dashboard_id`), length(`dashboard_id`) == 1)
        self$`dashboard_id` <- `dashboard_id`
      }
      if (!missing(`look_id`)) {
        stopifnot(is.numeric(`look_id`), length(`look_id`) == 1)
        self$`look_id` <- `look_id`
      }
      if (!missing(`content_type`)) {
        stopifnot(is.character(`content_type`), length(`content_type`) == 1)
        self$`content_type` <- `content_type`
      }
      if (!missing(`inherits`)) {
        self$`inherits` <- `inherits`
      }
      if (!missing(`inheriting_id`)) {
        stopifnot(is.numeric(`inheriting_id`), length(`inheriting_id`) == 1)
        self$`inheriting_id` <- `inheriting_id`
      }
      if (!missing(`can`)) {
        self$`can` <- `can`
      }
    },
    toJSON = function() {
      ContentMetaObject <- list()
      if (!is.null(self$`id`)) {
        ContentMetaObject[['id']] <- self$`id`
      }
      if (!is.null(self$`name`)) {
        ContentMetaObject[['name']] <- self$`name`
      }
      if (!is.null(self$`parent_id`)) {
        ContentMetaObject[['parent_id']] <- self$`parent_id`
      }
      if (!is.null(self$`space_id`)) {
        ContentMetaObject[['space_id']] <- self$`space_id`
      }
      if (!is.null(self$`dashboard_id`)) {
        ContentMetaObject[['dashboard_id']] <- self$`dashboard_id`
      }
      if (!is.null(self$`look_id`)) {
        ContentMetaObject[['look_id']] <- self$`look_id`
      }
      if (!is.null(self$`content_type`)) {
        ContentMetaObject[['content_type']] <- self$`content_type`
      }
      if (!is.null(self$`inherits`)) {
        ContentMetaObject[['inherits']] <- self$`inherits`
      }
      if (!is.null(self$`inheriting_id`)) {
        ContentMetaObject[['inheriting_id']] <- self$`inheriting_id`
      }
      if (!is.null(self$`can`)) {
        ContentMetaObject[['can']] <- self$`can`
      }

      ContentMetaObject
    },
    fromJSON = function(ContentMetaJson) {
      ContentMetaObject <- jsonlite::fromJSON(ContentMetaJson)
      if (!is.null(ContentMetaObject$`id`)) {
        self$`id` <- ContentMetaObject$`id`
      }
      if (!is.null(ContentMetaObject$`name`)) {
        self$`name` <- ContentMetaObject$`name`
      }
      if (!is.null(ContentMetaObject$`parent_id`)) {
        self$`parent_id` <- ContentMetaObject$`parent_id`
      }
      if (!is.null(ContentMetaObject$`space_id`)) {
        self$`space_id` <- ContentMetaObject$`space_id`
      }
      if (!is.null(ContentMetaObject$`dashboard_id`)) {
        self$`dashboard_id` <- ContentMetaObject$`dashboard_id`
      }
      if (!is.null(ContentMetaObject$`look_id`)) {
        self$`look_id` <- ContentMetaObject$`look_id`
      }
      if (!is.null(ContentMetaObject$`content_type`)) {
        self$`content_type` <- ContentMetaObject$`content_type`
      }
      if (!is.null(ContentMetaObject$`inherits`)) {
        self$`inherits` <- ContentMetaObject$`inherits`
      }
      if (!is.null(ContentMetaObject$`inheriting_id`)) {
        self$`inheriting_id` <- ContentMetaObject$`inheriting_id`
      }
      if (!is.null(ContentMetaObject$`can`)) {
        self$`can` <- ContentMetaObject$`can`
      }
    },
    toJSONString = function() {
       sprintf(
        '{
           "id": %d,
           "name": %s,
           "parent_id": %d,
           "space_id": %s,
           "dashboard_id": %s,
           "look_id": %d,
           "content_type": %s,
           "inherits": %s,
           "inheriting_id": %d,
           "can": %s
        }',
        self$`id`,
        self$`name`,
        self$`parent_id`,
        self$`space_id`,
        self$`dashboard_id`,
        self$`look_id`,
        self$`content_type`,
        self$`inherits`,
        self$`inheriting_id`,
        self$`can`
      )
    },
    fromJSONString = function(ContentMetaJson) {
      ContentMetaObject <- jsonlite::fromJSON(ContentMetaJson)
      self$`id` <- ContentMetaObject$`id`
      self$`name` <- ContentMetaObject$`name`
      self$`parent_id` <- ContentMetaObject$`parent_id`
      self$`space_id` <- ContentMetaObject$`space_id`
      self$`dashboard_id` <- ContentMetaObject$`dashboard_id`
      self$`look_id` <- ContentMetaObject$`look_id`
      self$`content_type` <- ContentMetaObject$`content_type`
      self$`inherits` <- ContentMetaObject$`inherits`
      self$`inheriting_id` <- ContentMetaObject$`inheriting_id`
      self$`can` <- ContentMetaObject$`can`
    }
  )
)

ext: r

max_stars_repo_name: adamlouis/lookr

max_issues_repo_name: adamlouis/lookr

max_forks_repo_name: adamlouis/lookr


text:
require( tidyverse )
require( viridis )
require( Seurat )
require( ComplexHeatmap)

# performs normalization, quality filtering of TaqMan assay cells.
# also builds heatmaps for all cells and for ltk positive cells only
#
# written by Seva Makeev, 2017 - 2021

plotDPI <- 600

workDir <- getwd()

resDir		 <- file.path( workDir, "Res")
dir.create( resDir, showWarnings = FALSE)

scTablesDir	<- file.path( resDir, "scTables")
dir.create( scTablesDir, showWarnings = FALSE)

taqmanSourceDir	 	<- file.path( workDir, "SourceData", "Taqman")

taqmanResDir	<- file.path( resDir, "taqman")
dir.create( taqmanResDir, showWarnings = FALSE)

plotDir 	<- file.path( resDir, "Plots")
dir.create( plotDir, showWarnings = FALSE)

taqmanPlotDir	<- file.path( plotDir, "taqman")
dir.create( taqmanPlotDir, showWarnings = FALSE)

taqmanPlotDir100dpi	<- file.path( taqmanPlotDir, "100dpi")
dir.create( taqmanPlotDir100dpi, showWarnings = FALSE)

taqmanPlotDir600dpi	<- file.path( taqmanPlotDir, "600dpi")
dir.create( taqmanPlotDir600dpi, showWarnings = FALSE)

plotHeight 	<- 10
plotWidth 	<- 10
Margin		<- 2

source("R/correctGeneNames.r")

taqmanData		<- read_delim( file = file.path( taqmanSourceDir, "30hpf_TAQMAN sum_paper.csv"), na = c("NA"), delim = ";")
colnames(taqmanData)	<- correctGeneNames( colnames(taqmanData))

taqmanGenes		<- setdiff( colnames( taqmanData), c("Sample", "rpl13", "Kan "))

taqmanData		<- as_tibble(sapply( taqmanData, function(x) ifelse( x %in% c("Undetermined", "undetermined", "Underemined", "underemined"), 40, x)))
taqmanData		<- type_convert( taqmanData)
taqmanData		<- rename( taqmanData, Kan = 'Kan ')

taqmanData 		<- mutate_at(taqmanData, setdiff( colnames( taqmanData), "Sample"), .funs = list( ~(40 - . )*log10(2)))
taqmanData		<- filter( taqmanData, Kan > quantile( apply(taqmanData[ , c("Kan")], 1, sum), 0.05))# re
taqmanData_Norm		<- mutate_at(taqmanData, .vars = vars(taqmanGenes), .funs = list( ~ .-Kan)) #we work in the log domain
#taqmanData_Norm		<- mutate_at(taqmanData, .vars = vars(taqmanGenes), .funs = list( ~ .-0.5*(Kan + rpl13))) #we work in the log domain
#taqmanData_Norm	<- taqmanData

taqmanMatrix	<- as.matrix(taqmanData_Norm[ , taqmanGenes])
rownames(taqmanMatrix) <- taqmanData$Sample

taqmanMatrix	<- apply( taqmanMatrix, 1:2, function(x) ifelse( is.na(x), 0, x))

keep <- apply(taqmanMatrix, 1, sum) > quantile( apply(taqmanMatrix, 1, sum), 0.05)
taqmanMatrixF <- taqmanMatrix[keep, ]

taqmanMatrixF <- taqmanMatrixF - min( taqmanMatrixF)

dens			<- density(as.matrix( taqmanMatrixF))
expLogMinimal		<- optimize(approxfun(dens$x,dens$y),interval=c(2,8))$minimum
taqmanMatrixF 		<- apply( taqmanMatrixF, 1:2, function(x) ifelse( x > expLogMinimal-1, x, 0))

densPlot	<- ggplot( data = data.frame( dens$x, dens$y), mapping = aes( dens$x, dens$y)) +
			geom_line() +
			geom_vline( xintercept = expLogMinimal - 1 , color = "red") +
			scale_x_continuous( name = "log10 probe expression" ) +
			scale_y_continuous( name = "Probe expression density") +
			theme(	plot.title 	= element_text( size = 50, hjust = 0),
				axis.title	= element_text( size = 30),
				axis.text	= element_text( size = 30),
				axis.text.y	= element_text( margin = margin( l = 0, t = 0, r = 1, b = 0))
				)
			#ggtitle( "taqman probes")

if (plotDPI == 600) {

  ggsave( paste0( "taqmanDensityPlot", ".png"), path = taqmanPlotDir600dpi, device = "png" , plot = densPlot, width = plotWidth, height = plotHeight, units = "cm", dpi = 600, scale = 4)

}

if (plotDPI == 100) {

  ggsave( paste0( "taqmanDensityPlot", ".png"), path = taqmanPlotDir600dpi, device = "png" , plot = densPlot, width = plotWidth, height = plotHeight, units = "cm", dpi = 100, scale = 4)

}

#now normalize for the minimum - which means substract the minimal value in the log scale

taqmanScTableDir 	<- file.path( scTablesDir, "taqman")
dir.create( taqmanScTableDir , showWarnings = FALSE)

write.table( taqmanMatrixF, file = file.path( taqmanScTableDir, "taqmanLogExps.csv"), sep = "\t")

source("R/seuratNorm.r")
source("R/getFinalClusterTypes.r")

clusterTreeDir		<- file.path( plotDir, "clusterTrees")
dir.create( clusterTreeDir, showWarnings = FALSE)

seuratTaqman 		<- seuratNorm( "taqman" )
genePos			<- sapply( taqmanGenes, function(x) WhichCells(seuratTaqman, subset.name = x, accept.low = 0))

Pigment			<- union( union( genePos[["ltk"]], genePos[["mitfa"]]), genePos[["pax7b"]])
phox2bbPigment		<- intersect( genePos[["phox2bb"]], Pigment)
phox2bbPigmentRat	<- length( phox2bbPigment)/length( genePos[["phox2bb"]])

neurog1Pigment		<- intersect( genePos[["neurog1"]], Pigment)
neurog1PigmentRat	<- length( neurog1Pigment)/length( genePos[["neurog1"]])

clusterOrderedCells	<- colnames( seuratTaqman@data[ , order( t(seuratTaqman@data)[ ,"sox10"], t(seuratTaqman@data)[ ,"mitfa"], t(seuratTaqman@data)[, "neurog1"], t(seuratTaqman@data)[, "phox2bb"], t(seuratTaqman@data)[, "ltk"], t(seuratTaqman@data)[, "pax7b"], decreasing = TRUE )] )
clusterOrderedCells	<- intersect( clusterOrderedCells, genePos["ltk"])
	#numbers are easier to sort; we will need this for the heatmap
levels( seuratTaqman@ident)	<- names( getFinalClusterTypes( seuratTaqman))

heatMapHeight 	<- 7
heatMapWidth 	<- 10
Margin		<- 2

source("R/drawTaqmanHeatMap.r")

if (plotDPI == 600) {
png( file = file.path( taqmanPlotDir600dpi, "taqmanBiclustHeatMap.png"),
	height = heatMapHeight + Margin + 1,
	width =  heatMapWidth + 2*Margin + 1,
	units = "in",
	res = plotDPI,
	pointsize = 2
)
	draw( drawTaqmanHeatMap( seuratTaqman,  clusterOrderedCells = colnames( seuratTaqman@data),  heatMapHeight, heatMapWidth, showCellNames = FALSE))
dev.off()}

if (plotDPI == 100) {
png( file = file.path( taqmanPlotDir100dpi, "taqmanBiclustHeatMap.png"),
	height = heatMapHeight + Margin + 1,
	width =  heatMapWidth + 2*Margin + 1,
	units = "in",
	res = plotDPI,
	pointsize = 2
)
	draw( drawTaqmanHeatMap( seuratTaqman,  clusterOrderedCells = colnames( seuratTaqman@data), heatMapHeight, heatMapWidth, showCellNames = FALSE))
dev.off()}

seuratTaqmanLtk <- SubsetData( seuratTaqman, cells.use = genePos$ltk)

if (plotDPI == 600) {
png( file = file.path( taqmanPlotDir600dpi, "taqmanLtkPosHeatMap.png"),
	height = heatMapHeight + Margin + 1,
	width =  heatMapWidth + 2*Margin + 1,
	units = "in",
	res = plotDPI,
	pointsize = 2
)
	draw( drawTaqmanHeatMap( seuratTaqmanLtk,  clusterOrderedCells = colnames( seuratTaqmanLtk@data),  heatMapHeight, heatMapWidth, showCellNames = FALSE))
dev.off()}

if (plotDPI == 100) {
png( file = file.path( taqmanPlotDir100dpi, "taqmanLtkPosHeatMap.png"),
	height = heatMapHeight + Margin + 1,
	width =  heatMapWidth + 2*Margin + 1,
	units = "in",
	res = plotDPI,
	pointsize = 2
)

	draw( drawTaqmanHeatMap( seuratTaqmanLtk,  clusterOrderedCells = colnames( seuratTaqmanLtk@data),  heatMapHeight, heatMapWidth, showCellNames = FALSE))

dev.off()}

png( filename = file.path( taqmanPlotDir, "taqmanGeneCors.png"))
	heatmap( cor(taqmanMatrixF), col = viridis(1024))
dev.off()

max_stars_repo_name: SevaVigg/NanostringDanioNCCscAnalysis

max_issues_repo_name: SevaVigg/NanostringDanioNCCscAnalysis

max_forks_repo_name: SevaVigg/NanostringDanioNCCscAnalysis


text:
"
a<- с(7:4, 0)
b<-с(8, 10.5, 0, -2, 9)
Написать программу, выполняющую следующие действия над векторами:
1) Сложение векторов a + b
2) Умножение векторов a * b
3) Деление векторов a / b
4) Нахождение среднего арифметического для каждого вектора
5) Нахождение суммы элементов каждого вектора
6) Представление результатов расчетов на экране
Указание. В программах использовать функции print(), paste(), paste0(), sum(), mean().
"

a <- c(7:4, 0)
print("Вектор А")
print(a)

b <- c(8, 10.5, 0, -2, 9)
print("Вектор B")
print(b)

result1 <- a + b
print(result1)

result2 <- a * b
print(result2)

result3 <- a / b
print(result3)

#Сред арифметическое
print(paste("Сред арифметическое А =",mean(a)))
print(paste("Сред арифметическое B =",mean(b)))

#Сумма элементов каждого вектора
print(paste0("Сумма элементов А=", sum(a)))
print(paste0("Сумма элементов B=", sum(b)))

max_stars_repo_path: Course II/R/pract/pract4/task2.r

max_stars_repo_name: GeorgiyDemo/FA

max_issues_repo_path: Course II/R/pract/pract4/task2.r

max_issues_repo_name: GeorgiyDemo/FA

max_forks_repo_path: Course II/R/pract/pract4/task2.r

max_forks_repo_name: GeorgiyDemo/FA


text:
# /* !/usr/bin/Rscript */

#' ## Model selection
#' Compare and assess best-of-breed models to determine candidate with
#' optimal fit for predicting survival of the Titanic disaster.
#' Best-of-breed is decided by initial output metrics, e.g. lowest
#' deviance score between GLM models, or Accuracy for other types.
#'
#' ### Receiver Operating Characteristic Curves
#' ROC curves plot the *Sensitivity* (true positive rate) of a model
#' against the *Specificity* (false positive rate) at various threshold
#' settings. In conjunction with the summary statistic AUC (Area Under
#' Curve) it can be used to assess the performance of a model.
#'
#' The following snippet plots ROC curves for a number of models
#' trained earlier in order to directly compare their performance
#' against one another. This is accomplished using the `plotROC`
#' function from the *discern* package which wraps the process of
#' creating confusion matrices to generate the plots and arranging
#' them on one canvas. The package can be installed from:
#' `install_github('andybeeching/discern')`
plotROC(
  models    = list(glmFit2, glmFit9, rfFit, cfFit, adaFit),
  res       = "Fate",
  cv        = cv,
  colors    = c("orange", "green", "blue", "yellow", "red"),
  labels    = c("glmFit2", "glmFit9", "rfFit", "cfFit", "adaFit")
)

#' ### Re-sampling Result Accuracy Distribution
#' Another way to assess models is to compare the distribution of
#' results from their respective k-fold validations. Repeated
#' validation means each model has 30 results - one for every
#' re-sample.
#'
#' To compare, gather the results for each model through caret's
#' `resample` method, then plot.
cvResults <- resamples(list(
  glm1 = glmFit1,
  glm2 = glmFit2,
  glm8 = glmFit8,
  glm9 = glmFit9,
  ada  = adaFit,
  rf   = rfFit,
  rf2  = rfFit,
  cf   = cfFit,
  svm  = svmFit
))
summary(cvResults)
dotplot(cvResults, metric="ROC")
bwplot(cvResults) # from lattice

#' ### Learning curves
#' An important performance facet to examine for any model is its overall
#' generalisability, that is, determine whether it may be inherently biased
#' (overfit), or conversely contain too much variance (underfit).
#'
#' One way to assert bias vs variance in a model is to plot how well
#' it performs when trained on increasingly larger amounts of data.
#'
#' Essentially the training dataset is broken into sub-training sets
#' based on ratio intervals (e.g. 10%, 20%, 30% etc...), and each
#' time the model is trained (with the original formula) it is used
#' to predict the response variable of the known cross-validation
#' dataset, and the Accuracy recorded. The resultant line is known as
#' a "learning curve".
#'
#' Ideally the curve would oscillate a little at the start as the
#' amount of data used to train the model doubles and triples, but
#' then flatten out as observations normalise and fill-out the
#' observation space.
#'
#' Once more a method from the "discern" package is used to perform
#' the above analysis. `plotLearningCurves` iteratively trains a list
#' of models using the cIn the below example on GLM models are
plotLearningCurves(
  models = list(glmFit2, glmFit9),
  labels = c("Glm2", "Glm9"),
  metric = "ROC",
  ctrlFn = fitControl,
  df     = tr,
  cv     = cv,
  colors = c("orange", "blue"),
  seed   = SEED
)

# /*
writeLines("\n-------------")
writeLines("Analysis: DONE")
# */

max_stars_repo_name: andybeeching/kaggle-titanic

max_issues_repo_name: andybeeching/kaggle-titanic

max_forks_repo_name: andybeeching/kaggle-titanic


text:
# Before running this code, you should install following packages:

#devtools::install_github("rstudio/sparklyr")
#install.packages('broom')

library(sparklyr)

config <- spark_config()
config[["spark.r.command"]] <- "./r_env.zip/r_env/bin/Rscript"
config[["spark.yarn.dist.archives"]] <- "r_env.zip"
config$sparklyr.apply.env.R_HOME <- "./r_env.zip/r_env/lib/R"
config$sparklyr.apply.env.RHOME <- "./r_env.zip/r_env"
config$sparklyr.apply.env.R_SHARE_DIR <- "./r_env.zip/r_env/lib/R/share"
config$sparklyr.apply.env.R_INCLUDE_DIR <- "./r_env.zip/r_env/lib/R/include"

sc <- spark_connect(master = "yarn-client", config = config)

sdf_len(sc, 5, repartition = 1) %>%
  spark_apply(function(e) I(e))

iris_tbl <- sdf_copy_to(sc, iris)

spark_apply(
  iris_tbl,
  function(e) broom::tidy(lm(Petal_Length ~ Petal_Width, e)),
  names = c("term", "estimate", "std.error", "statistic", "p.value"),
  group_by = "Species",
  packages = FALSE)

max_stars_repo_name: chezou/sparklyr-distribute

max_issues_repo_name: chezou/sparklyr-distribute

max_forks_repo_name: chezou/sparklyr-distribute


text:
### Preparation of QHI phenology data for the Coastal Phenology Analysis
## 8 Feb 2018 Jakob Assmann j.assmann@ed.ac.uk

# dependencies
library(dplyr)

# Housekeeping
script_path <- "phenology_data/qhi/"

qhi_phen <- read.csv("phenology_data/qhi/qiki_phen_with_before_2017.csv")

### The aim is to create a df with the following structure and
### only containing: ERIVAG flowering, SALARC green_up and DRYINT flowering
###
# * site_name (SITE)
# * spp (SPECIES)
# * phen_stage (flowering, green_up)
# * site_spp_phase (SITE_SPECIES_phen)
# * spp_phase (SPECIES_phen)
# * plot_id (plot_id)
# * year (YYYY)
# * doy (J)
# * prior_obs (J)
# * snowmelt (J)

### ERIVAG flowering
erivag_flowering <- qhi_phen %>% filter(Spp == "ERIVAG") %>%
  select (Spp, Plot.ID, Year, P2, P2_before, P1)
erivag_flowering$phen_stage <- "flowering"
names(erivag_flowering) <- c("spp", "plot_id", "year", "doy", "prior_visit", "snowmelt", "phen_stage")
erivag_flowering$site_name <- "QHI"
erivag_flowering <- erivag_flowering %>%
  mutate(spp_phen = paste0(spp, "_", phen_stage), site_spp_phen = paste0(site_name, "_", spp, "_", phen_stage) )
erivag_flowering <- erivag_flowering[, c(8,1,7,10,9,2,3,4,5,6)]

### SALARC green_up
salarc_green_up <- qhi_phen %>% filter(Spp == "SALARC") %>%
  select (Spp, Plot.ID, Year, P2, P2_before, P1)
salarc_green_up$phen_stage <- "green_up"
names(salarc_green_up) <- c("spp", "plot_id", "year", "doy", "prior_visit", "snowmelt", "phen_stage")
salarc_green_up$site_name <- "QHI"
salarc_green_up <- salarc_green_up %>%
  mutate(spp_phen = paste0(spp, "_", phen_stage), site_spp_phen = paste0(site_name, "_", spp, "_", phen_stage) )
salarc_green_up <- salarc_green_up[, c(8,1,7,10,9,2,3,4,5,6)]

### DRYINT flowering
dryint_flowering <- qhi_phen %>% filter(Spp == "DRYINT") %>%
  select (Spp, Plot.ID, Year, P2, P2_before, P1)
dryint_flowering$phen_stage <- "flowering"
names(dryint_flowering) <- c("spp", "plot_id", "year", "doy", "prior_visit", "snowmelt", "phen_stage")
dryint_flowering$site_name <- "QHI"
dryint_flowering <- dryint_flowering %>%
  mutate(spp_phen = paste0(spp, "_", phen_stage), site_spp_phen = paste0(site_name, "_", spp, "_", phen_stage) )
dryint_flowering <- dryint_flowering[, c(8,1,7,10,9,2,3,4,5,6)]

### Check consistency
names(erivag_flowering) == names(salarc_green_up)
names(salarc_green_up) == names(dryint_flowering)

### Combine into one dataframe
qhi_phen <- rbind(erivag_flowering, salarc_green_up, dryint_flowering)

### Create unique plot identifiers
qhi_phen$plot_id <- paste0(qhi_phen$spp, "_", qhi_phen$plot_id)

### Check whether there are intevals wider than 10
sum(qhi_phen$doy - qhi_phen$prior_visit > 10, na.rm = T)
qhi_phen[which(qhi_phen$doy - qhi_phen$prior_visit > 10),]
mean(qhi_phen[which(qhi_phen$doy - qhi_phen$prior_visit > 10),]$doy -
       qhi_phen[which(qhi_phen$doy - qhi_phen$prior_visit > 10),]$prior_visit)
# 30 day mean window for these two data points, ach well...
# wonder what happened that year?

### Export data
write.csv(qhi_phen, paste0(script_path, "qhi_phen.csv"), row.names=FALSE)

max_stars_repo_name: jakobjassmann/coastalphenology

max_issues_repo_name: jakobjassmann/coastalphenology

max_forks_repo_name: jakobjassmann/coastalphenology


text:
cat(sapply(strsplit(readLines(tail(commandArgs(), n=1)), " "), function(s) {
  paste(Reduce(function(x,y) {
    if (x[1] == 0) {
      c(1, y)
    } else {
      if (x[length(x)] == y) {
        x[length(x)-1] <- as.integer(x[length(x)-1]) + 1
        x
      } else {
        append(x, c(1, y))
      }
    }
  }, s, c(0, "")), collapse=" ")
}), sep="\n")

max_stars_repo_name: IlkhamGaysin/ce-challenges

max_issues_repo_name: IlkhamGaysin/ce-challenges

max_forks_repo_name: IlkhamGaysin/ce-challenges


text:
#
# Example R code to install packages
# See http://cran.r-project.org/doc/manuals/R-admin.html#Installing-packages for details
#

install.packages("e1071", dependencies = TRUE)


text:
#' Extract genotype data from PLINK .raw files for Illumina 450K SNPs
#'
#' @param filenames A vector of filenames of PLINK .raw files from which to extract
#' genotype data.
#' @return Matrix with rows corresponding to SNPs, columns to samples and values
#' equal to 0, 1 or 2 corresponding to genotypes.
#' @examples
#' R> writeLines(meffil.snp.names("450k"), con="snp-names.txt")
#' shell> plink --bfile dataset --extract snp-names.txt --recodeA --out genotypes.raw --noweb
#' R> filenames <- "genotypes.raw"
#' R> genotypes <- meffil.extract.genotypes(filenames)
#'
#' @export
meffil.extract.genotypes <- function(filenames, verbose=F) {
    stopifnot(all(sapply(filenames, file.exists)))

    table.list <- lapply(filenames, function(filename) {
        msg("Reading plink file", filename, verbose=verbose)
        read.table(filename, header=T)
    })

    sample.names <- table.list[[1]][,1] ## family id
    genotypes <- lapply(table.list, function(genotype.table) {
        genotype.table[match(sample.names, genotype.table[,1]), ## match family ids
                       -(1:6),
                       drop=F]
    })
    genotypes <- do.call(cbind, genotypes)
    colnames(genotypes) <- sub("_.*", "", colnames(genotypes))
    rownames(genotypes) <- sample.names

    t(as.matrix(genotypes))
}

max_stars_repo_name: RichardJActon/meffil

max_issues_repo_name: RichardJActon/meffil

max_forks_repo_name: RichardJActon/meffil


text:
1
-- @@stderr --
dtrace: invalid probe specifier 0: invalid probe identifier 0: Unknown probe identifier

max_stars_repo_name: alan-maguire/dtrace-utils

max_issues_repo_name: tjfontaine/dtrace-utils

max_forks_repo_name: tjfontaine/dtrace-utils


text:
################################################################################

GbmTrain <- function(train.name="train",target=NULL,args=list(test.name=NULL,
	R=1,pls=FALSE,prune=FALSE,std=FALSE,features=c("static","alsfrs","fvc",
	"svc","vital"), slope.features=NULL)) {
	# Trains a GBM model given training data.
	# Args:
	#	train.name: base name of training dataset
	#	target: target vector; if NULL, future slopes are used as target
	#	args: list of additional arguments
	#	-test.name: base name of test set, if you want training set to be
	#	 processed using both test and train data
	#	- pls: TRUE means predicting the last score rather than the slope
	#	- prune: TRUE means keeping only features with data for at least 918/2
	#		patients
	#	- std: TRUE means standardize columns by subtracting away column means
	#	- features: a character vector of feature sets to get for the training
	#		data.  Choose from "static","alsfrs","fvc","svc","vital","lab"
	#	- slope.features: a character vector of datasets to get slope features
	#		for.  Choose from NULL if no slopes are desired, or among "static",
	#		"alsfrs","fvc","svc","vital","lab"
	#
	# Returns:
	#	a function that given a test set, returns predictions
	#
	# Side effects:
	#	Saves predictions to file
	#
	# Example usage:
	if(FALSE) {
		features = c("static","alsfrs","fvc","svc","vital")
	   	slope.features = NULL
	   	args=list(test.name="test",R=1,pls=FALSE,prune=FALSE,std=FALSE,
	   		features=features,slope.features=slope.features)
		model <- GbmTrain(train.name="train",target=NULL,args=args)
		predictions<-model("test")
	}

	# Preprocess training data according to arguments
	train.preprocess<-PreprocessFeatures(train.name=train.name,
		test.name=args[["test.name"]],args=args)

	# Get training data and target
	train.data<-as.data.frame(train.preprocess$feature.data)

	# Need to throw out if not imputing NAs
	#train.data<-train.data[,colnames(train.data)!="slope.height.slope"]

	target<-train.preprocess$target

	# Fit multiple models to training data
	GetModel <- function(i) {
		# Fit to training data
		object <- gbm(target~.,data=train.data[,-1],distribution="gaussian",n.trees=20000,
			cv.folds=10,verbose=FALSE)

		# Show features used
		print(summary(object,plot=FALSE))

		# Find optimal number of iterations
		#best.iter<-gbm.perf(object,method="OOB") # Check with OOB
		#best.iter<-gbm.perf(object,method="test") # Check with 50% heldout test
		best.iter<-gbm.perf(object,method="cv",plot.it=FALSE) # Check with CV

		# Report summary results
		print(sprintf("Model: %d, CV error: %g, Best iter: %g",i,min(object$cv.error),best.iter))
		return(list(object = object, param = best.iter))
	}
	models <- lapply(1:args[["R"]], GetModel)

	GbmPredict <- function(test.name="test",save.preds=TRUE) {
		# Returns predictions based on the trained model given a test set
		#
		# Args:
		#	test.name: base name of test dataset
		#   save.preds: save predictions to file? TRUE by default
		#
		# Returns:
		#	Two column matrix, ordered by the first column, with first column
		#   "subject.id" containing test set subject ids and the second column
		#   "prediction" containing predicted slopes.
		#
		# Side effects:
		#	Saves predictions to file

		# Preprocess test data according to arguments
		test.preprocess<-PreprocessFeatures(train.name=NULL,test.name=test.name,
			args=args,train.cols=train.preprocess$train.cols,train.data.std=
			train.preprocess$train.data.std,
			train.data.na=train.preprocess$train.data.na)

		test.data<-as.data.frame(test.preprocess$feature.data)

		# Need to throw out if not imputing NAs
		#test.data<-test.data[,colnames(test.data)!="slope.height.slope"]

		predictions <- rowMeans(sapply(models, function(model)
			predict(model$object,test.data[,-1],model$param)))

		# If predicting last score, transform predictions back to slopes
		if (args[["pls"]]) {
			scaling<-PlsScale(test.name)
			predictions<-predictions-
				scaling$score.first*scaling$time.scale
		}
		predictions <- cbind(test.data[,1],predictions)

		if(save.preds) {
			# Write predictions to file
			out.file = GetPredictionFilename("GbmTrain",train.name,test.name,args)
			WritePredictions(predictions, out.file)
		}
		return(predictions)
	}
	return(GbmPredict)
}


text:
#Caching
#Joe Swintek June 13,2013
#######################################################################################
#										On Lists				     				  #
#######################################################################################
rm(list = ls()) #Clears the workspace
Cache<-list()   #Declare list
Cache

Cache[[1]] #out of bounds error

Cache$A #Will give a null

is.null(Cache$A) #Will give a TRUE

Cache$A<-c(10,13) #Give Cache A some Values
Cache$A
Cache[[1]]  #Usefull oddity

Cache$Z<-'This is Cache Z' #Give Cache Z a string value
Cache$B<-list('Empty'={}) #Give Cache B a value of an empty list
Cache
#Will associate indexes with the order at which it was assigned
Cache[[2]]  #Same as Cache Z
Cache[[3]]  #Same as Cache B

length(Cache) #Number of top most elements ,i.e. 3

Cache[[4]]  #Gives an out of bounds error

#Treating a list as a Hash Table
Cache[['A']] #Same as Cache$A
Cache[['Z']] #Same as Cache$Z
Cache[['B']] #Same as Cache$B

Key<-'Z'

Cache[[Key]] #will print Cache Z
Cache[['B']][['Empty']]<-'The label is a misnomer'
Cache
Cache[['B']]

#######################################################################################
#									Recursive Fibonacci				     			  #
#######################################################################################

RecursiveFibonacci<-function(n){ #Slow
	#This gives the Fibonacci number after n steps
    if (n<=2) return(1)
return (RecursiveFibonacci(n-1)+RecursiveFibonacci(n-2))
}

RecursiveFibonacci(1)   #Should be 1
RecursiveFibonacci(2)   #Should be 1
RecursiveFibonacci(7)   #Should be 13
RecursiveFibonacci(20)  #Should be 6765
RecursiveFibonacci(30)  #takes some time

system.time(RecursiveFibonacci(30)) #Time it takes
system.time(RecursiveFibonacci(32))
#######################################################################################
#						Aliasing functions (Well... kind of)					   	  #
#######################################################################################

PrintFunction<-function(Function,Args){
	#This will run a function with arguments Args
	print(Function(Args))
return()
}

PrintFunction(RecursiveFibonacci,10) #prints 55
Function<-RecursiveFibonacci
Function(10) #prints 55

#######################################################################################
#									On Caching Itself				     			  #
#######################################################################################

CacheFunction<-function(Cache,FunctionAsString,Args){
#This function will draw a function result from Cache or run the function and Cache the result
#This will return the updated Cache and print the result
#FunctionAsString must be string
#Cache is the Cache being used
#Args is a vector of arguments for the function

	Key<-FunctionAsString  #Gets the Cache Key
	Function<-get(FunctionAsString)  #Turn a sting into code of its name
									 #You can use get to write a function that writes a function

	#Check to see if that Function is Cached with the args
	InCache<-!is.null(Cache[[Key]][[toString(Args)]]) #is.null checks for a null argument
													  #toString will convert a vector to a string of its elements
													  #The ! means not in R

	#If the function/Args is not in the Cache run it and Cache the function
	if (InCache==FALSE){
		#if (!is.null(Cache[[Key]])==FALSE){Cache[[Key]]$CasheSetUp<-TRUE} #Fixes index problem
		Cache[[Key]][[toString(Args)]]<-Function(Args)
	}

	#Print the Result
	print(Cache[[Key]][[toString(Args)]])

return (Cache) #Return the Cache
}

Cache<-list()  #Reset Cache
system.time(Cache<-CacheFunction(Cache,'RecursiveFibonacci',30)) #Run the function and store the Cache
system.time(Cache<-CacheFunction(Cache,'RecursiveFibonacci',30)) #Rerun and it works!

#Test
system.time(Cache<-CacheFunction(Cache,'RecursiveFibonacci',32)) #Fails to Run...
print(Cache)

#Retry
Cache<-list()  #Reset Cache
system.time(Cache<-CacheFunction(Cache,'RecursiveFibonacci',30)) #Run the function and store the Cache
system.time(Cache<-CacheFunction(Cache,'RecursiveFibonacci',32)) #Works this time

print(Cache)

#######################################################################################
#										Global Caching 						   		  #
#######################################################################################
CacheFunction<-function(FunctionAsString,Args){
#This function will draw a function result from Cache or run the function and Cache the result
#This will return the updated Cache and print the result
#FunctionAsString must be string
#Cache is the Cache being used And is a global
#Args is a vector of arguments for the function

	Key<-FunctionAsString  #Gets the Cache Key
	Function<-get(FunctionAsString)  #Turn a sting into code of its name
									 #You can use get to write a function that writes a function

	#Check to see if that Function is Cached with the args
	InCache<-!is.null(Cache[[Key]][[toString(Args)]]) #is.null checks for a null argument
													  #toString will convert a vector to a string of its elements
													  #The ! means not in R
	#If the function/Args is not in the Cache run it and Cache the function
	if (InCache==FALSE){
		if (!is.null(Cache[[Key]])==FALSE){Cache[[Key]]$CasheSetUp<<-TRUE} #Fixes index problem
		Cache[[Key]][[toString(Args)]]<<-Function(Args) #<<-Assigns the variable to the upper most workspace
	}

	#Print the Result
return(Cache[[Key]][[toString(Args)]])
}

Cache<-list()  #Reset Cache
system.time(Number<-CacheFunction('RecursiveFibonacci',30)) #Run the function and store the Cache
print(Number)

system.time(Number<-CacheFunction('RecursiveFibonacci',30)) #Rerun and it works!
print(Number)

print(Cache)

system.time(Number<-CacheFunction('RecursiveFibonacci',33)) #Try it with 33
print(Number)

system.time(Number<-CacheFunction('RecursiveFibonacci',33)) #Rerun and it works!
print(Number)

print(Cache)

#######################################################################################
#								Save and load the Cache						   		  #
#######################################################################################

save('Cache',file="MyCache.RData") #Saves Cache to file
ls()						#look at workspace
rm('Cache')					#Delete the Cache
ls()
load("MyCache.RData")  		#loads the Cashe into the workspace
					        #This will overide any variables that have the same name as ones in the loaded file
ls()  					    #look at workspace to see what loaded

CacheSize<-object.size(Cache) #lists the size of the data object
print(CacheSize,units='auto') #Prints the Cache size in different units

#######################################################################################
#								Final Demonstration 								  #
#######################################################################################

Cache<-list()  #Reset Cache
RecursiveFibonacci<-function(n){ #Slow, but now cached for speed
	#This gives the Fibonacci number after n steps
    if (n<=2) return(1)
return (CacheFunction('RecursiveFibonacci',(n-1))+CacheFunction('RecursiveFibonacci',(n-2)))
}

system.time(print(RecursiveFibonacci(30))) #Time it takes
system.time(print(RecursiveFibonacci(33)))
system.time(print(RecursiveFibonacci(100)))
system.time(print(RecursiveFibonacci(200)))

length(Cache)

CacheSize<-object.size(Cache) #lists the size of the data object
print(CacheSize,units='auto') #Prints the Cache size in different units

max_stars_repo_name: hbrittingham/R-User-Group

max_issues_repo_name: hbrittingham/R-User-Group

max_forks_repo_name: hbrittingham/R-User-Group


text:
getwd()
setwd("/mnt/workspace/rstudio_projects/data/")

library(dplyr)

load('gdmt_bb.RData')
view(df_gdmt_bb)
names(df_gdmt_bb)

#labels----------------------------------------------------
df_gdmt_bb$ACEIdisc <- factor(df_gdmt_bb$ACEIdisc,
                            levels = c(1,2,3),
                            labels = c('Yes- prescribed', 'No- not prescribed', 'No- contraindicated'))

df_gdmt_bb$AHA_ARB <- factor(df_gdmt_bb$AHA_ARB,
                           levels = c(2,3,4,5,6,7,8,9),
                           labels = c('candesartan', 'irbesartan', 'losartan', 'valsartan', 'telmisartan', 'eprosartan',
                                       'tasosartan', 'other'))

df_gdmt_bb$ARBdisc <- factor(df_gdmt_bb$ARBdisc,
                           levels = c(1,2,3),
                           labels = c('Yes- prescribed', 'No- not prescribed', 'No- contraindicated'))

df_gdmt_bb$JC_ACESPE <- factor(df_gdmt_bb$JC_ACESPE,
                             levels = c(2,4,5,9,26,29,32,35,47,50,51),
                             labels = c('quinapril', 'perindopril', 'ramipril', 'captopril', 'benazepril', 'trandolapril',
                                        'fosinopril', 'lisinopril', 'enalopril', 'lisinopril', 'other'))

df_gdmt_bb$aha_arni <- factor(df_gdmt_bb$aha_arni,
                            levels = c(1,2,3),
                            labels = c('sacubitril/valsartan', 'other', 'unknown'))

#subset out vars of interest----------------------------------------------------
df_raasi <- subset(df_gdmt_bb, select = c(PATIENT_ID,
                                          ACEIdisc, JC_ACEDISC, JC_ACESPE,
                                          ACEDISCCONTRASPEC_Missing, AHA_ACEDISCCONTRASPEC_1:AHA_ACEDISCCONTRASPEC_5, JC_ACEDISCCONTRA,
                                          AHA_ARB, ARBdisc, JC_ARBDISC,
                                          AHA_ARBDISCCONTRASPEC_1:AHA_ARBDISCCONTRASPEC_5, ARBDISCCONTRASPEC_Missing, JC_ARBDISCCONTRA,
                                          aha_arni, jc_arnidisc,
                                          aha_arnidisccontraspec_1:aha_arnidisccontraspec_8, jc_arnidisccontra,
                                          aha_arniaceiarb,
                                          HF_ACEIARBforLVSD))

#create variables summarizing contraindications----------------------------------------------------

#ACEI contraindicated? binary var
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ACEIcontra = ifelse(AHA_ACEDISCCONTRASPEC_1 == 1 | AHA_ACEDISCCONTRASPEC_2 == 1 | AHA_ACEDISCCONTRASPEC_3 == 1|
                                  AHA_ACEDISCCONTRASPEC_4 == 1 | AHA_ACEDISCCONTRASPEC_5 == 1, 1, 0), .after=AHA_ACEDISCCONTRASPEC_5)
df_raasi$GDMT_ACEIcontra[is.na(df_raasi$GDMT_ACEIcontra)] <- 0
addmargins(table(df_raasi$GDMT_ACEIcontra, exclude = NULL))

#ACEI contraindicated? -> expanded to include those patients with contra as designated by ACEIdisc
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ACEIcontra_exp = ifelse(GDMT_ACEIcontra == 0 & ACEIdisc == 'No- contraindicated', 1, GDMT_ACEIcontra), .after=GDMT_ACEIcontra)
df_raasi$GDMT_ACEIcontra_exp[is.na(df_raasi$GDMT_ACEIcontra_exp)] <- 0
addmargins(table(df_raasi$GDMT_ACEIcontra, df_raasi$GDMT_ACEIcontra_exp, exclude = NULL))

#ARB contraindicated? binary var
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARBcontra = ifelse(AHA_ARBDISCCONTRASPEC_1 == 1 | AHA_ARBDISCCONTRASPEC_2 == 1 | AHA_ARBDISCCONTRASPEC_3 == 1 |
                                 AHA_ARBDISCCONTRASPEC_4 == 1 | AHA_ARBDISCCONTRASPEC_5 == 1, 1, 0), .after=AHA_ARBDISCCONTRASPEC_5)
df_raasi$GDMT_ARBcontra[is.na(df_raasi$GDMT_ARBcontra)] <- 0
addmargins(table(df_raasi$GDMT_ARBcontra, exclude = NULL))

#ACEI contraindicated? -> expanded to include those patients with contra as designated by ACEIdisc
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARBcontra_exp = ifelse(GDMT_ARBcontra == 0 & ARBdisc == 'No- contraindicated', 1, GDMT_ARBcontra), .after=GDMT_ARBcontra)
#df_raasi$GDMT_ARBcontra_exp[is.na(df_raasi$GDMT_ARBcontra_exp)] <- 0
addmargins(table(df_raasi$GDMT_ARBcontra, df_raasi$GDMT_ARBcontra_exp, exclude = NULL))
addmargins(table(df_raasi$ARBdisc, df_raasi$GDMT_ARBcontra_exp, exclude = NULL))

#ARNI contraindicated? binary var
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARNIcontra = ifelse(aha_arnidisccontraspec_1 == 1 | aha_arnidisccontraspec_2 == 1 | aha_arnidisccontraspec_3 == 1 |
                                  aha_arnidisccontraspec_4 == 1 | aha_arnidisccontraspec_5 == 1 | aha_arnidisccontraspec_6 == 1 |
                                  aha_arnidisccontraspec_7 == 1 | aha_arnidisccontraspec_8 == 1, 1, 0), .after=aha_arnidisccontraspec_8)
df_raasi$GDMT_ARNIcontra[is.na(df_raasi$GDMT_ARNIcontra)] <- 0
addmargins(table(df_raasi$GDMT_ARNIcontra, exclude = NULL))

#tables----------------------------------------------------
names(df_raasi)
lapply(subset(df_raasi, select = -c(PATIENT_ID, AHA_ACEDISCCONTRASPEC_1:AHA_ACEDISCCONTRASPEC_5, ACEDISCCONTRASPEC_Missing,
                                    AHA_ARBDISCCONTRASPEC_1:AHA_ARBDISCCONTRASPEC_5, ARBDISCCONTRASPEC_Missing,
                                    aha_arnidisccontraspec_1:aha_arnidisccontraspec_8)),
                                 function(x) addmargins(table(x, exclude = NULL)))

#ACEI#####################################################
#part I: comparing the contra-indication variables
addmargins(table(df_raasi$ACEIdisc, df_raasi$GDMT_ACEIcontra, exclude = NULL))
  #2 elements here: 1) 538pts for whom ACEI is contra got it anyways
  #                 2) disturbingly there are 7036 patients for whom ACEI based on ACEdisc was contra-indicated who did not
  #                    register in the ACEIcontra variable > may need to expand it
addmargins(table(df_raasi$ACEIdisc, df_raasi$GDMT_ACEIcontra_exp, exclude = NULL))
  #here's the follow up comparison of the expansion var - the expanded version is valid
addmargins(table(df_raasi$ACEDISCCONTRASPEC_Missing, df_raasi$GDMT_ACEIcontra, exclude = NULL))
  #valid - all patients with a missing contraindication do not have a subcategory contra listed

#part II: who actually got an ACEI on dc?
addmargins(table(df_raasi$ACEIdisc, df_raasi$JC_ACEDISC, exclude = NULL))
  #ACEIdisc and JC_ACEDISC are consistent but ACEIdisc is better, has less missing and more info - differentiates no/contra from no
addmargins(table(df_raasi$JC_ACESPE, df_raasi$ACEIdisc, exclude = NULL))

#ARB#####################################################
#part III: comparing the contra-indication variables
addmargins(table(df_raasi$ARBdisc, df_raasi$GDMT_ARBcontra, exclude = NULL))
  #512 patients had a contraindicaiton to ARB and still got one
  #7230 patients without a condraindication designated by AHA_ARBDISCCONTRASPEC_ vars actually had a contraindication per ARBdisc; build an exp form?
addmargins(table(df_raasi$ARBdisc, df_raasi$GDMT_ARBcontra_exp, exclude = NULL))
  #follow up has folded the 7230 pts noted above into the contraindicated pts for the GDMT_ARBcontra var
addmargins(table(df_raasi$ARBDISCCONTRASPEC_Missing, df_raasi$GDMT_ARBcontra, exclude = NULL))
  #it seems as aif the _Missing variable is derived from the HA_ARBDISCCONTRASPEC_ var series
addmargins(table(df_raasi$ARBdisc, df_raasi$JC_ARBDISCCONTRA, exclude = NULL))

#part IV: who's actually getting an ARB on dc?
addmargins(table(df_raasi$ARBdisc, df_raasi$AHA_ARB, exclude = NULL))
  #there are 2 patients who received losartan who are not counted as having received an ARB, otherwise ARBdisc is a very good variable
addmargins(table(df_raasi$ARBdisc, df_raasi$JC_ARBDISC, exclude = NULL))
  #JC_ARBDISC is a good variable but there are some some values designated by NA that ARBdisc has knowledge of

#ARNI#####################################################
#part V: comparing the contra-indications for ARNI
addmargins(table(df_raasi$jc_arnidisc, df_raasi$GDMT_ARNIcontra, exclude = NULL))
addmargins(table(df_raasi$aha_arni, df_raasi$GDMT_ARNIcontra, exclude = NULL))
  #85+8 patients with contra to ARNI got one anyways

addmargins(table(df_raasi$GDMT_ARNIcontra, df_raasi$jc_arnidisccontra, exclude = NULL))
addmargins(table(df_raasi$jc_arnidisc, df_raasi$jc_arnidisccontra, exclude = NULL))

#part VI: who actually got an ARNI?
addmargins(table(df_raasi$jc_arnidisc, exclude = NULL))
addmargins(table(df_raasi$aha_arni, exclude = NULL))
addmargins(table(df_raasi$jc_arnidisc, df_raasi$aha_arni, exclude = NULL))
  #8679 pts got an ARNI that the aha_arni classification didn't know about??
  #1+3 pts got an ARNI that jc_arni classificaiotn didn't know about - we will give the coders the benefit of the doubt

#other checks
#who's getting both ACEI and ARB??
addmargins(table(df_raasi$ARBdisc, df_raasi$ACEIdisc, exclude = NULL))
  #there are 772 patient's prescribed both??
addmargins(table(df_raasi$AHA_ARB, df_raasi$JC_ACESPE, exclude = NULL))
  #there are 1+5 + 2 +1 +16 73 19 + 1 + 2 + 9+ 2 + 18  patient's getting an ACEI/ARB combo that can be identified??
#who's getting ARNI/ACEI or ARNI/ARB?
addmargins(table(df_raasi$ARBdisc, df_raasi$aha_arni, exclude = NULL))
  #there are 526 patients prescribed ARNI and ARB?
addmargins(table(df_raasi$ACEIdisc, df_raasi$aha_arni, exclude = NULL))
  #there are 183 patients prescribed ACEI and ARNI on dc?

#create essential variables summarizing prescribed and prescribed/contra----------------------------------------------------

#did the patient get an ACEI on dc?
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ACEI = ifelse(df_raasi$JC_ACESPE %in% 'lisinopril' | df_raasi$ACEIdisc == 'Yes- prescribed', 1,
                     ifelse(df_raasi$ACEIdisc == 'No- not prescribed' | df_raasi$ACEIdisc == 'No- contraindicated', 0, 99)))
addmargins(table(df_raasi$GDMT_ACEI, df_raasi$ACEIdisc, exclude = NULL))

#this test case of %in% vs == is very interesting... I'm still not sure why R is doing what it's doing
#df_raasi <-
#  df_raasi %>%
#  mutate(GDMT_ACEI = ifelse(df_raasi$JC_ACESPE %in% 'lisinopril' | df_raasi$ACEIdisc %in% 'Yes- prescribed', 1,
#                            ifelse(df_raasi$ACEIdisc %in% 'No- not prescribed' | df_raasi$ACEIdisc %in% 'No- contraindicated', 0, 99)))

#did the patient get an ACEI on dc or was it contraindicated?
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ACEIorContra = ifelse(df_raasi$GDMT_ACEI == 1 | df_raasi$GDMT_ACEIcontra_exp == 1, 1, 0))
addmargins(table(df_raasi$GDMT_ACEIorContra, df_raasi$ACEIdisc, exclude = NULL))

#did the patient get an ARB on dc? # 34822 + 2
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARB = ifelse(df_raasi$ARBdisc == 'Yes- prescribed' | df_raasi$AHA_ARB %in% 'losartan', 1,
                    ifelse(df_raasi$ARBdisc == 'No- not prescribed' | df_raasi$ARBdisc == 'No- contraindicated', 0, 99)))
addmargins(table(df_raasi$GDMT_ARB, df_raasi$ARBdisc, exclude = NULL))

#did the patient get an ARB on dc or was it contraindicated?
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARBorContra = ifelse(df_raasi$GDMT_ARB == 1 | df_raasi$GDMT_ARBcontra_exp == 1, 1, 0))
addmargins(table(df_raasi$GDMT_ARBorContra, df_raasi$ARBdisc, exclude = NULL))

#did the patient get an ARNI on dc? #22657
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARNI = ifelse(df_raasi$jc_arnidisc == 1 | df_raasi$aha_arni %in% c('sacubitril/valsartan', 'other'), 1, 0))
addmargins(table(df_raasi$GDMT_ARNI, exclude = NULL))

#did the patient get an ARNI on dc or was it contraindicated?
df_raasi <-
  df_raasi %>%
  mutate(GDMT_ARNIorContra = ifelse(df_raasi$GDMT_ARNI == 1 | df_raasi$GDMT_ARNIcontra == 1, 1, 0))
addmargins(table(df_raasi$GDMT_ARNIorContra, exclude = NULL))

#implement the key var into the gdmt dataframe----------------------------------------------------

#ACEI contraindicated? binary var
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ACEIcontra = ifelse(AHA_ACEDISCCONTRASPEC_1 == 1 | AHA_ACEDISCCONTRASPEC_2 == 1 | AHA_ACEDISCCONTRASPEC_3 == 1|
                                  AHA_ACEDISCCONTRASPEC_4 == 1 | AHA_ACEDISCCONTRASPEC_5 == 1, 1, 0))
df_gdmt_bb$GDMT_ACEIcontra[is.na(df_gdmt_bb$GDMT_ACEIcontra)] <- 0
addmargins(table(df_gdmt_bb$GDMT_ACEIcontra, exclude = NULL))

#ACEI contraindicated? -> expanded to include those patients with contra as designated by ACEIdisc
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ACEIcontra_exp = ifelse(GDMT_ACEIcontra == 0 & ACEIdisc == 'No- contraindicated', 1, GDMT_ACEIcontra))
df_gdmt_bb$GDMT_ACEIcontra_exp[is.na(df_gdmt_bb$GDMT_ACEIcontra_exp)] <- 0
addmargins(table(df_gdmt_bb$GDMT_ACEIcontra, df_gdmt_bb$GDMT_ACEIcontra_exp, exclude = NULL))

#ARB contraindicated? binary var
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARBcontra = ifelse(AHA_ARBDISCCONTRASPEC_1 == 1 | AHA_ARBDISCCONTRASPEC_2 == 1 | AHA_ARBDISCCONTRASPEC_3 == 1 |
                                 AHA_ARBDISCCONTRASPEC_4 == 1 | AHA_ARBDISCCONTRASPEC_5 == 1, 1, 0))
df_gdmt_bb$GDMT_ARBcontra[is.na(df_gdmt_bb$GDMT_ARBcontra)] <- 0
addmargins(table(df_gdmt_bb$GDMT_ARBcontra, exclude = NULL))

#ACEI contraindicated? -> expanded to include those patients with contra as designated by ACEIdisc
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARBcontra_exp = ifelse(GDMT_ARBcontra == 0 & ARBdisc == 'No- contraindicated', 1, GDMT_ARBcontra))
#df_raasi$GDMT_ARBcontra_exp[is.na(df_raasi$GDMT_ARBcontra_exp)] <- 0
addmargins(table(df_gdmt_bb$GDMT_ARBcontra, df_gdmt_bb$GDMT_ARBcontra_exp, exclude = NULL))
addmargins(table(df_gdmt_bb$ARBdisc, df_gdmt_bb$GDMT_ARBcontra_exp, exclude = NULL))

#ARNI contraindicated? binary var
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARNIcontra = ifelse(aha_arnidisccontraspec_1 == 1 | aha_arnidisccontraspec_2 == 1 | aha_arnidisccontraspec_3 == 1 |
                                  aha_arnidisccontraspec_4 == 1 | aha_arnidisccontraspec_5 == 1 | aha_arnidisccontraspec_6 == 1 |
                                  aha_arnidisccontraspec_7 == 1 | aha_arnidisccontraspec_8 == 1, 1, 0))
df_gdmt_bb$GDMT_ARNIcontra[is.na(df_gdmt_bb$GDMT_ARNIcontra)] <- 0
addmargins(table(df_gdmt_bb$GDMT_ARNIcontra, exclude = NULL))

#did the patient get an ACEI on dc?
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ACEI = ifelse(JC_ACESPE %in% 'lisinopril' | ACEIdisc == 'Yes- prescribed', 1,
                     ifelse(ACEIdisc == 'No- not prescribed' | ACEIdisc == 'No- contraindicated', 0, 99)))
addmargins(table(df_gdmt_bb$GDMT_ACEI, df_gdmt_bb$ACEIdisc, exclude = NULL))

#did the patient get an ACEI on dc or was it contraindicated?
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ACEIorContra = ifelse(GDMT_ACEI == 1 | GDMT_ACEIcontra_exp == 1, 1, 0))
addmargins(table(df_gdmt_bb$GDMT_ACEIorContra, df_gdmt_bb$ACEIdisc, exclude = NULL))

#did the patient get an ARB on dc? # 34822 + 2
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARB = ifelse(ARBdisc == 'Yes- prescribed' | AHA_ARB %in% 'losartan', 1,
                    ifelse(ARBdisc == 'No- not prescribed' | ARBdisc == 'No- contraindicated', 0, 99)))
addmargins(table(df_gdmt_bb$GDMT_ARB, df_gdmt_bb$ARBdisc, exclude = NULL))

#did the patient get an ARB on dc or was it contraindicated?
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARBorContra = ifelse(GDMT_ARB == 1 | GDMT_ARBcontra_exp == 1, 1, 0))
addmargins(table(df_gdmt_bb$GDMT_ARBorContra, df_gdmt_bb$ARBdisc, exclude = NULL))

#did the patient get an ARNI on dc? #22657
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARNI = ifelse(jc_arnidisc == 1 | aha_arni %in% c('sacubitril/valsartan', 'other'), 1, 0))
addmargins(table(df_gdmt_bb$GDMT_ARNI, exclude = NULL))

#did the patient get an ARNI on dc or was it contraindicated?
df_gdmt_bb <-
  df_gdmt_bb %>%
  mutate(GDMT_ARNIorContra = ifelse(GDMT_ARNI == 1 | GDMT_ARNIcontra == 1, 1, 0))
addmargins(table(df_gdmt_bb$GDMT_ARNIorContra, exclude = NULL))

#clean up excess dc Rx BB-related var----------------------------------------------------
names(df_gdmt_bb)
df_gdmt_bb_raasi <- subset(df_gdmt_bb, select = -c(ACEIdisc, JC_ACEDISC, JC_ACEDISCCONTRA, ACEDISCCONTRASPEC_Missing,
                                                   JC_ARBDISC, ARBDISCCONTRASPEC_Missing, JC_ARBDISCCONTRA,
                                                   aha_arni, jc_arnidisc, jc_arnidisccontra,
                                                   aha_arniaceiarb,
                                                   HF_ACEIARBforLVSD))
names(df_gdmt_bb_raasi)

#save the new data frame----------------------------------------------------
#save(df_gdmt_bb_raasi, file = '/mnt/workspace/rstudio_projects/gdmt_bb_raasi.RData')


text:
{{rimport}}('plot.r', '__init__.r', 'sampleinfo.r')

exprfile  = {{i.expr      | R}}
batchfile = {{i.batch     | R}}
tool      = {{args.tool    | R}}
outfile   = {{o.outfile  | R}}
outdir    = {{o.outdir   | R}}
plot      = {{args.plot    | R}}
ggs       = {{args.ggs     | R}}
devpars   = {{args.devpars | R}}
hmrows    = {{args.hmrows  | R}}

expr    = read.table.nodup (exprfile, sep = "\t", header = TRUE, row.names = 1, check.names = F)

samples    = colnames(expr)
sampleinfo = SampleInfo(batchfile)
batch      = sampleinfo$dataframe(sampleinfo$select(sample = samples, get = c('Sample', 'Batch')))

batch      = factor(batch[samples, 'Batch'])

if (tool == 'combat') {
	library(sva)
	newexpr   = ComBat(dat=expr, batch=batch, par.prior = TRUE, mod = NULL)
	write.table (round(newexpr, 3), outfile, col.names=T, row.names=T, sep="\t", quote=F)
} else { # leave it for further tools
	stop('Unsupported tool: {{args.tool}}.')
}

# boxplot
if (plot$boxplot) {
	bpfile = file.path(outdir, "{{i.expr | fn2}}.boxplot.png")
	plot.boxplot(newexpr, bpfile, stack = T, devpars = devpars, ggs = ggs$boxplot)
}

# heatmap
if (plot$heatmap) {
	hmfile = file.path(outdir, "{{i.expr | fn2}}.heatmap.png")
	hmexp  = if (nrow(newexpr) > hmrows) newexpr[sample(nrow(newexpr),size=hmrows),] else newexpr
	plot.heatmap(hmexp, hmfile, devpars = devpars, ggs = ggs$heatmap)
}

# histgram
if (plot$histogram) {
	histfile = file.path(outdir, "{{i.expr | fn2}}.histo.png")
	plot.histo(as.data.frame(newexpr), histfile, devpars = devpars, ggs = ggs$histogram)
}

max_stars_repo_name: LeaveYeah/bioprocs

max_issues_repo_name: LeaveYeah/bioprocs

max_forks_repo_name: LeaveYeah/bioprocs


text:
# Looker API 3.0 Reference
#
# ### Authorization  The Looker API uses Looker **API3** credentials for authorization and access control. Looker admins can create API3 credentials on Looker's **Admin/Users** page. Pass API3 credentials to the **/login** endpoint to obtain a temporary access_token. Include that access_token in the Authorization header of Looker API requests. For details, see [Looker API Authorization](https://looker.com/docs/r/api/authorization)  ### Client SDKs  The Looker API is a RESTful system that should be usable by any programming language capable of making HTTPS requests. Client SDKs for a variety of programming languages can be generated from the Looker API's Swagger JSON metadata to streamline use of the Looker API in your applications. A client SDK for Ruby is available as an example. For more information, see [Looker API Client SDKs](https://looker.com/docs/r/api/client_sdks)  ### Try It Out!  The 'api-docs' page served by the Looker instance includes 'Try It Out!' buttons for each API method. After logging in with API3 credentials, you can use the \"Try It Out!\" buttons to call the API directly from the documentation page to interactively explore API features and responses.  ### Versioning  Future releases of Looker will expand this API release-by-release to securely expose more and more of the core power of Looker to API client applications. API endpoints marked as \"beta\" may receive breaking changes without warning. Stable (non-beta) API endpoints should not receive breaking changes in future releases. For more information, see [Looker API Versioning](https://looker.com/docs/r/api/versioning)
#
# OpenAPI spec version: 3.0.0
#
# Generated by: https://github.com/swagger-api/swagger-codegen.git

#' ContentFavorite Class
#'
#' @field id
#' @field user_id
#' @field content_metadata_id
#' @field look_id
#' @field dashboard_id
#' @field look
#' @field dashboard
#' @field can
#'
#' @importFrom R6 R6Class
#' @importFrom jsonlite fromJSON toJSON
#' @export
ContentFavorite <- R6::R6Class(
  'ContentFavorite',
  public = list(
    `id` = NULL,
    `user_id` = NULL,
    `content_metadata_id` = NULL,
    `look_id` = NULL,
    `dashboard_id` = NULL,
    `look` = NULL,
    `dashboard` = NULL,
    `can` = NULL,
    initialize = function(`id`, `user_id`, `content_metadata_id`, `look_id`, `dashboard_id`, `look`, `dashboard`, `can`){
      if (!missing(`id`)) {
        stopifnot(is.numeric(`id`), length(`id`) == 1)
        self$`id` <- `id`
      }
      if (!missing(`user_id`)) {
        stopifnot(is.numeric(`user_id`), length(`user_id`) == 1)
        self$`user_id` <- `user_id`
      }
      if (!missing(`content_metadata_id`)) {
        stopifnot(is.numeric(`content_metadata_id`), length(`content_metadata_id`) == 1)
        self$`content_metadata_id` <- `content_metadata_id`
      }
      if (!missing(`look_id`)) {
        stopifnot(is.numeric(`look_id`), length(`look_id`) == 1)
        self$`look_id` <- `look_id`
      }
      if (!missing(`dashboard_id`)) {
        stopifnot(is.numeric(`dashboard_id`), length(`dashboard_id`) == 1)
        self$`dashboard_id` <- `dashboard_id`
      }
      if (!missing(`look`)) {
        stopifnot(R6::is.R6(`look`))
        self$`look` <- `look`
      }
      if (!missing(`dashboard`)) {
        stopifnot(R6::is.R6(`dashboard`))
        self$`dashboard` <- `dashboard`
      }
      if (!missing(`can`)) {
        self$`can` <- `can`
      }
    },
    toJSON = function() {
      ContentFavoriteObject <- list()
      if (!is.null(self$`id`)) {
        ContentFavoriteObject[['id']] <- self$`id`
      }
      if (!is.null(self$`user_id`)) {
        ContentFavoriteObject[['user_id']] <- self$`user_id`
      }
      if (!is.null(self$`content_metadata_id`)) {
        ContentFavoriteObject[['content_metadata_id']] <- self$`content_metadata_id`
      }
      if (!is.null(self$`look_id`)) {
        ContentFavoriteObject[['look_id']] <- self$`look_id`
      }
      if (!is.null(self$`dashboard_id`)) {
        ContentFavoriteObject[['dashboard_id']] <- self$`dashboard_id`
      }
      if (!is.null(self$`look`)) {
        ContentFavoriteObject[['look']] <- self$`look`$toJSON()
      }
      if (!is.null(self$`dashboard`)) {
        ContentFavoriteObject[['dashboard']] <- self$`dashboard`$toJSON()
      }
      if (!is.null(self$`can`)) {
        ContentFavoriteObject[['can']] <- self$`can`
      }

      ContentFavoriteObject
    },
    fromJSON = function(ContentFavoriteJson) {
      ContentFavoriteObject <- jsonlite::fromJSON(ContentFavoriteJson)
      if (!is.null(ContentFavoriteObject$`id`)) {
        self$`id` <- ContentFavoriteObject$`id`
      }
      if (!is.null(ContentFavoriteObject$`user_id`)) {
        self$`user_id` <- ContentFavoriteObject$`user_id`
      }
      if (!is.null(ContentFavoriteObject$`content_metadata_id`)) {
        self$`content_metadata_id` <- ContentFavoriteObject$`content_metadata_id`
      }
      if (!is.null(ContentFavoriteObject$`look_id`)) {
        self$`look_id` <- ContentFavoriteObject$`look_id`
      }
      if (!is.null(ContentFavoriteObject$`dashboard_id`)) {
        self$`dashboard_id` <- ContentFavoriteObject$`dashboard_id`
      }
      if (!is.null(ContentFavoriteObject$`look`)) {
        lookObject <- LookBasic$new()
        lookObject$fromJSON(jsonlite::toJSON(ContentFavoriteObject$look, auto_unbox = TRUE))
        self$`look` <- lookObject
      }
      if (!is.null(ContentFavoriteObject$`dashboard`)) {
        dashboardObject <- DashboardBase$new()
        dashboardObject$fromJSON(jsonlite::toJSON(ContentFavoriteObject$dashboard, auto_unbox = TRUE))
        self$`dashboard` <- dashboardObject
      }
      if (!is.null(ContentFavoriteObject$`can`)) {
        self$`can` <- ContentFavoriteObject$`can`
      }
    },
    toJSONString = function() {
       sprintf(
        '{
           "id": %d,
           "user_id": %d,
           "content_metadata_id": %d,
           "look_id": %d,
           "dashboard_id": %d,
           "look": %s,
           "dashboard": %s,
           "can": %s
        }',
        self$`id`,
        self$`user_id`,
        self$`content_metadata_id`,
        self$`look_id`,
        self$`dashboard_id`,
        self$`look`$toJSON(),
        self$`dashboard`$toJSON(),
        self$`can`
      )
    },
    fromJSONString = function(ContentFavoriteJson) {
      ContentFavoriteObject <- jsonlite::fromJSON(ContentFavoriteJson)
      self$`id` <- ContentFavoriteObject$`id`
      self$`user_id` <- ContentFavoriteObject$`user_id`
      self$`content_metadata_id` <- ContentFavoriteObject$`content_metadata_id`
      self$`look_id` <- ContentFavoriteObject$`look_id`
      self$`dashboard_id` <- ContentFavoriteObject$`dashboard_id`
      LookBasicObject -> LookBasic$new()
      self$`look` <- LookBasicObject$fromJSON(jsonlite::toJSON(ContentFavoriteObject$look, auto_unbox = TRUE))
      DashboardBaseObject -> DashboardBase$new()
      self$`dashboard` <- DashboardBaseObject$fromJSON(jsonlite::toJSON(ContentFavoriteObject$dashboard, auto_unbox = TRUE))
      self$`can` <- ContentFavoriteObject$`can`
    }
  )
)


text:
library(ggplot2)
library(lme4)
library(stringr)
library(dplyr)
library(tidyr)
library(scales)
library(readr)
library(tidyverse)
library(fs)
library(purrr)
library(svglite)

blog_width = 2400
blog_height =1000
blog_dpi = 175

windowsFonts("OpenSans" = windowsFont("OpenSans"))
aligners = c('fave', 'maus', 'mfa_default_train', 'mfa_ipat_train', 'mfa_enlish_adapt', 'mfa_enlish', 'mfa_english_ipa', 'mfa_english_ipa_adapt')

directory = 'D:\\Data\\speech\\benchmark_datasets\\buckeye\\accuracy_data'

cbbPalette <- c("#FB5607", "#8338EC", "#FF006E", "#FFBE0B",  "#3A86FF", '#6EC200', '#BA3200')
theme_memcauliffe <- function () {
  base_size = 12
  half_line <- base_size/2
  theme_bw(base_size=base_size, base_family="OpenSans") %+replace%
    theme(
      plot.background = element_rect(fill = "#003566"),
      panel.background = element_rect(fill = "#003566", color = '#003566'),
      panel.border = element_rect(color='#000814', fill='NA'),
      panel.grid.major = element_line(size = 0.25, linetype = 'solid',
                                      colour = "#FFC300"),
      panel.grid.minor = element_line(size = 0.125, linetype = 'solid',
                                      colour = "#FFD60A"),
      axis.text = element_text(size = rel(1.2), colour = "#FFD60A"),
      axis.text.y = element_text(margin = margin(r = 0.8*half_line/2),
                                 hjust = 1, colour = "#FFD60A"),
      axis.title = element_text(colour = "#FFC300", size = rel(1.6)),
      plot.title = element_text(colour = "#FFC300", size = rel(1.8), hjust=0,
                                margin = margin(b = half_line * 1.2)),
      strip.background = element_rect(fill = "#FFC300", colour = NA),
      strip.text = element_text(colour = "#000814", size = rel(1.2),  margin = margin( b = 3, t = 3,  l = 4, r = 4)),
      legend.background = element_blank(),
      legend.key = element_blank(),
      legend.title = element_text(size = rel(1.2), colour = "#FFC300"),
      legend.text = element_text(size = rel(1), colour = "#FFC300"),
    )
}

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  require(dplyr)

  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }

  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  dots <- lapply(groupvars, as.symbol)

  functions <- list("n()", paste("mean(",measurevar,")", sep = ''), as.formula(paste("~median(",measurevar,")", sep = '')), paste("sd(",measurevar,")", sep = ''))

  # hack to deal with median specified as formula
  names <-  lapply(functions, function(x){if(length(x)>1){as.character(x[2])}else{x[1]}})
  names <- lapply(names, as.name)

  datac <- data %>% group_by_(.dots = dots) %>% summarise_(.dots = functions) %>% rename_(.dots = setNames(names, c('N', 'mean', 'median', 'sd')))

  # Rename the "mean" column

  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval:
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult

  return(datac)
}

### Phonsay

#phonsay <- read.csv('results.txt')
#phonsay <- subset(phonsay, prec_consonant %in% c('B','D','G','F','S','P','T','K','HH'))
#phonsay <- subset(phonsay, foll_consonant %in% c('B','D','G','F','S','P','T','K','HH'))

#phonsay <- subset(phonsay, word != 'carpet')
#phonsay$vowel <- factor(str_replace(phonsay$vowel, '[0-2]',''))

#phonsay$cons_type <- 'Sonorant'
#phonsay[phonsay$foll_consonant %in% c('D','Z','B','G'),]$cons_type <- 'Voiced'
#phonsay[phonsay$foll_consonant %in% c('T','S','P','K'),]$cons_type <- 'Voiceless'
#phonsay$cons_type <- factor(phonsay$cons_type, levels = c('Sonorant','Voiced','Voiceless'))

#phonsay$aligner <- factor(phonsay$aligner, labels = c('FAVE','MFA', 'PLA'))
#phonsay$data <- factor(phonsay$data, labels = c('Flat', 'Lab', 'LS', 'LS-Clean','SCOTUS'))

#phonsay$condition <- factor(paste(phonsay$aligner,phonsay$data,sep='-\n'))
#phonsay$diff_mag <- abs(phonsay$diff)

#phonsay$diff.std <- scale(phonsay$diff)

#phonsay$type <- factor(phonsay$type, levels = c('vowel_begin', 'vowel_end', 'cons_end'), labels = c('Vowel begin point', 'Vowel end point', 'Obstruent end point'))
#phonsay$condition <- factor(phonsay$condition, levels = c('mfa-flat', 'mfa-lab', 'mfa-librispeech', 'prosodylab-flat', 'prosodylab-lab'), labels = c('MFA-Flat', 'MFA-Lab', 'MFA-LibriSpeech', 'Prosodylab-Flat', 'Prosodylab-Lab'))

max_stars_repo_name: mmcauliffe/memcauliffe-blog-scripts

max_issues_repo_name: mmcauliffe/memcauliffe-blog-scripts

max_forks_repo_name: mmcauliffe/memcauliffe-blog-scripts


text:
/*
 ==============================================================================

 This file is part of the iPlug 2 library. Copyright (C) the iPlug 2 developers.

 See LICENSE.txt for  more info.

 ==============================================================================
*/

#include "config.h"   // This is your plugin's config.h.
#include <AudioUnit.r>

#define UseExtendedThingResource 1

#include <CoreServices/CoreServices.r>

// this is a define used to indicate that a component has no static data that would mean
// that no more than one instance could be open at a time - never been true for AUs
#ifndef cmpThreadSafeOnMac
#define cmpThreadSafeOnMac  0x10000000
#endif

#undef  TARGET_REZ_MAC_X86
#if defined(__i386__) || defined(i386_YES)
  #define TARGET_REZ_MAC_X86        1
#else
  #define TARGET_REZ_MAC_X86        0
#endif

#undef  TARGET_REZ_MAC_X86_64
#if defined(__x86_64__) || defined(x86_64_YES)
  #define TARGET_REZ_MAC_X86_64     1
#else
  #define TARGET_REZ_MAC_X86_64     0
#endif

#if TARGET_OS_MAC
   #if TARGET_REZ_MAC_X86 && TARGET_REZ_MAC_X86_64
    #define TARGET_REZ_FAT_COMPONENTS_2 1
    #define Target_PlatformType     platformIA32NativeEntryPoint
    #define Target_SecondPlatformType platformX86_64NativeEntryPoint
  #elif TARGET_REZ_MAC_X86
    #define Target_PlatformType     platformIA32NativeEntryPoint
  #elif TARGET_REZ_MAC_X86_64
    #define Target_PlatformType     platformX86_64NativeEntryPoint
  #else
    #error you gotta target something
  #endif
  #define Target_CodeResType    'dlle'
  #define TARGET_REZ_USE_DLLE   1
#else
  #error get a real platform type
#endif // not TARGET_OS_MAC

#ifndef TARGET_REZ_FAT_COMPONENTS_2
  #define TARGET_REZ_FAT_COMPONENTS_2   0
#endif

#ifndef TARGET_REZ_FAT_COMPONENTS_4
  #define TARGET_REZ_FAT_COMPONENTS_4   0
#endif

// ----------------

//#ifdef _DEBUG
//  #define PLUG_PUBLIC_NAME PLUG_NAME "_DEBUG"
//#else
#define PLUG_PUBLIC_NAME PLUG_NAME
//#endif

#define RES_ID 1000
#define RES_NAME PLUG_MFR ": " PLUG_PUBLIC_NAME

#if PLUG_TYPE==0
#if PLUG_DOES_MIDI_IN
#define COMP_TYPE kAudioUnitType_MusicEffect
#else
#define COMP_TYPE kAudioUnitType_Effect
#endif
#elif PLUG_TYPE==1
#define COMP_TYPE kAudioUnitType_MusicDevice
#elif PLUG_TYPE==2
#define COMP_TYPE 'aumi'
#endif

resource 'STR ' (RES_ID, purgeable) {
  RES_NAME
};

resource 'STR ' (RES_ID + 1, purgeable) {
  PLUG_PUBLIC_NAME " AU"
};

resource 'dlle' (RES_ID) {
  AUV2_ENTRY_STR
};

resource 'thng' (RES_ID, RES_NAME) {
  COMP_TYPE,
  PLUG_UNIQUE_ID,
  PLUG_MFR_ID,
  0, 0, 0, 0,               //  no 68K
  'STR ', RES_ID,
  'STR ', RES_ID + 1,
  0,  0,      // icon
  PLUG_VERSION_HEX,
  componentHasMultiplePlatforms | componentDoAutoVersion,
  0,
  {
    cmpThreadSafeOnMac,
    Target_CodeResType, RES_ID,
    Target_PlatformType,
#if TARGET_REZ_FAT_COMPONENTS_2 || TARGET_REZ_FAT_COMPONENTS_4
    cmpThreadSafeOnMac,
    Target_CodeResType, RES_ID,
    Target_SecondPlatformType,
#endif
  }
};

#undef RES_ID


text:
# -- 1/ create data :
	# We have 3 species : Dicoccoides / Dicoccum / Durum
	# We have 2 treatments : High or Low Nitrogen
	# We have 4 samples per specie per treatment.
	# For each sample we measure the expression of 5 genes
sample=paste(rep("sample_",24) , seq(1,24) , sep="")
specie=c(rep("dicoccoides" , 8) , rep("dicoccum" , 8) , rep("durum" , 8))
data=data.frame(sample,specie)
for (i in seq(1:5)){
  gene=c( sample(c(1:40),8) , sample(c(10:45),8) , sample(c(5:15),8))
  data=cbind(data , gene)
  colnames(data)[ncol(data)]=paste("gene_",i,sep="")
 }

#I arrange data to have a difference between low and High
data[data$specie=="durum" , c(3:7)]=data[data$specie=="durum" , c(3:7)]-30
data

# -- 2/ Calculate distances between samples
rownames(data)=data[,1]
distance=dist(data[ , c(3:7)] , diag=TRUE)

# -- 3/ Make the MDS
fit <- cmdscale(as.dist(distance) , eig=TRUE, k=2)
x <- fit$points[,1]
y <- fit$points[,2]

library(wordcloud)
png("#34_MDS_with_no_overlapping_names.png" , width = 480, height = 480)
par(mar=c(4,4,2,2))
textplot(x, y, data$sample, xlab="Coordinate 1", ylab="Coordinate 2", main="", col="white" , xlim=c(-65,78) )
#text(x,y, data$sample , col=as.numeric(data$specie) , lwd=5 , cex=1.3)
legend("topright", legend = c("durum","dicoccum","dicoccoides") , col = c("green","red","black") , pch =20 , bty = "n",  pt.cex = 1.7, cex = 1 , text.col = "black", horiz = FALSE, inset = c(0.1, 0.1))
dev.off()

max_stars_repo_name: JedStephens/R-graph-gallery

max_issues_repo_name: JedStephens/R-graph-gallery

max_forks_repo_name: JedStephens/R-graph-gallery


text:
N <- as.matrix(read.table("nTags.txt"))
###### input format for data
###### AB  A  B
###### AB: pet count between anchors A and B
###### A: tags in anchor A
###### B: tags in anchor B
data <- as.matrix(read.table("data.txt"))
nRows <- length(data[,1])
q <- data[,1] - 1
m <- data[,2]
n <- N - m
k <- data[,3]
x <- phyper(q, m, n, k, lower.tail = FALSE)

###### p-value adjustment with Benjamini-Hockberg method
fdr <- p.adjust(x, "BH")
a <- matrix(1, nRows, 4)
a[,1] <- x
a[,2] <- fdr
a[,3:4] <- -log10(a[,1:2]) ###### -log10(p-value)
a[is.infinite(a[,3]),3] <- 1000 ###### replace the extreme values with 1000
a[is.infinite(a[,4]),4] <- 1000 ###### replace the extreme values with 1000
a <- as.data.frame(a)
for (i in 1:2) {
  a[[i]] <- as.numeric(format(a[[i]], digits = 3))
}
for (i in 3:4) {
  a[[i]] <- round(a[[i]], 2)
}
write.table(a, file="result.txt", sep = "\t", row.names = FALSE, col.names = FALSE)


text:
library(magrittr)

# likelihoods of the three gems
a = 1 / 2
b = 1 / 3
c = 1 / 6

Ea = (1 +
      + b * (1 / c - 1) * a / (a + b)
      + c * (1 / b - 1) * a / (a + c)) / (1 - a)

Eb = (c * 1  + a * (1 + (1 / c - 1) * a / (a + b))) / (1 - b)
Ec = (b * 1  + a * (1 + (1 / b - 1) * a / (a + c))) / (1 - c)

a * Ea + b * Eb + c * Ec

#Simulate to verify.
gemCountSim <- function() {
  vecCount = c(0,0,0)
  while (TRUE) {
    draw = rmultinom(1,1,c(a,b,c)) %>% as.vector
    vecCount = vecCount + draw
    if (vecCount %>% min > 0)
      return(vecCount[1])
  }
}

replicate(10000,gemCountSim(),simplify = "array") %>% mean


text:
# load required libraries
library(tidyverse)
library(unmarked)
library(MASS)
library(parallel)
library(extraDistr)
library(binom)
library(abind)

# location of functions
source("functions.r")

# define input parameters
genYears <- 20 # number of years for a generation
demogV <- 0.4373 # demographic stochasticity variance - this is variance in individual annual growth rate among individuals - taken from Rhodes et al. (2011) based on the radio tacking data only
                # model. Here the standard deviation of the growth rate parameter estimate was 0.04458435 from 220 individuals, so we calculated demogV as (0.04458435^2) * 220
envtlV <- c(0, 0.01) # environmental stochasticity variance
envtlCorr <- c(0, 0.5, 1) # correlation in environmental stochasticity among populations
nSimsPop <- 1000 # number of replicates for the population simulations
nYears <- 50 # simulation time horizon
decline <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50) # percentage decline over three generations
f0 <- 0.0285 # distance detection parameter from Dique et al(2003), eff(W) = 1/f0
f0se <- 0.0021 # standard error of f0 from Dique et al(2003)
stripMissLow <- 0 # low end of range of proportion missed for strip transects from Dique et al. (2001)
stripMissMean <- 0.16 # mean proportion missed for strip transects from Dique et al. (2001)
stripMissHigh <- 0.33 # high end of range of proportion missed for strip transects from Dique et al. (2001)
tranDayLine <- 4600 # line transect length per day in m
tranDayArea <- 7.5 # area search area per day in ha
siteSize <- 100 # site size in ha
monitRep <- 21 # monitoring time frame
budget <- c(95, 190, 380) # field days per year
survIntens <- c(30, 60, 90) # proportion of the site surveyed
method <- c("line", "area") # "line" is all line and "area" is line transects outside the urban footprint and all area searches inside the urban footprint
monInter <- c(1, 2) # monitoring interval in years
strat <- c("equal", "area", "density", "invdensity")

# load strata and initial population densities - set -9999 to NA and remove strata with < 100 ha of habitat
# then group by genetic population and nest
# note here that we replace the density data for "core - urban footprint - north west" stratum with the density data for "core - urban footprint - north west" stratum
# since this density data is missing for "core - urban footprint - north west". we also use data from Biolink 2019 for North Stradbroke Island (NSI) and assume core and
# non-core habitat densities are the same for NSI. this is the file "start_dens_modified.csv"
pops <- read.csv("input/start_dens_modified.csv") %>%
  mutate(KD_MEAN = replace(KD_MEAN, KD_MEAN == -9999, NA), KD_SD = replace(KD_SD, KD_SD==-9999, NA),
    KDLN_MEAN = replace(KDLN_MEAN, KDLN_MEAN == -9999, NA), KDLN_SD=replace(KDLN_SD, KDLN_SD == -9999, NA)) %>%
  filter(AREA >= 100) %>% mutate(ABUND = round(KD_MEAN * AREA)) %>% as_tibble() %>% rowid_to_column("ID") %>% group_by(POP) %>% nest()
# calculate the population size for each genetic populations
pops <- pops %>% mutate(N = map_dbl(data, ~sum(.$ABUND)))

# simulate population dynamics - comment out if already run

procTime <- system.time({

popSimComb <- expand.grid(decline = decline, envtlV = envtlV, envtlCorr = envtlCorr)
# remove rows when environmental variation is 0 and correlation is 0 or 0.5
popSimComb <- popSimComb %>% filter(!(envtlV == 0 & ((envtlCorr == 0) | (envtlCorr == 0.5))))
cl <- makeCluster(detectCores() - 1)
clusterExport(cl, list("pops", "projPops", "nSimsPop", "nYears", "genYears", "demogV" ))
clusterEvalQ(cl, library(tidyverse))
clusterEvalQ(cl, library(MASS))
parApply(cl = cl, X = popSimComb, MARGIN = 1, FUN = function(x) {saveRDS(object = projPops(pops = pops$N, names = pops$POP, nYears = nYears, decline = x["decline"], genYears = genYears,
                demogV = demogV, envtlV = x["envtlV"], envtlCorr = x["envtlCorr"], nSims = nSimsPop), file = paste("output/popsims/popsim", "_decline",
                x["decline"], "_envtlV", x["envtlV"], "_envtlCorr", x["envtlCorr"], ".rds", sep = ""))})
stopCluster(cl)

})

print(procTime)

# generate survey designs - comment out if already run

procTime <- system.time({

# get combinations of parameters
survSimComb <- expand.grid(budget = budget, survIntens = survIntens, method = method, monInter = monInter, strat = strat, monitRep = monitRep)
# collate fixed survey parameters together
fixedSurvParams <- c(f0, f0se, tranDayLine, tranDayArea, siteSize)
names(fixedSurvParams) <- c("f0", "f0se", "tranDayLine", "tranDayArea", "siteSize")
cl <- makeCluster(detectCores() - 1)
clusterExport(cl, list("getSurveys", "smart.round", "pops", "fixedSurvParams" ))
clusterEvalQ(cl, library(tidyverse))
parApply(cl = cl, X = survSimComb, MARGIN = 1, FUN = function(x) {saveRDS(object = getSurveys(pops = pops, survParams = c(x, fixedSurvParams)), file = paste("output/surveys/survey", "_budget",
                x["budget"], "_survIntens", as.numeric(x["survIntens"]), "_method", x["method"], "_monInter", as.numeric(x["monInter"]), "_strat", x["strat"], "_monitRep", as.numeric(x["monitRep"]), ".rds", sep = ""))})
stopCluster(cl)

})

print(procTime)

# generate simulated survey data
# get survey information
Surveys <- apply(survSimComb, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/surveys/survey", "_budget", x["budget"], "_survIntens",
              as.numeric(x["survIntens"]), "_method", x["method"], "_monInter", as.numeric(x["monInter"]), "_strat", x["strat"], "_monitRep",
              as.numeric(x["monitRep"]), ".rds", sep = ""))})
# loop through population combinations
for (i in 1:nrow(popSimComb)) {

  procTime <- system.time({

  # get population simulations
  popSim <- readRDS(file = paste("output/popsims/popsim", "_decline", popSimComb[i, "decline"], "_envtlV",
        popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))
  cl <- makeCluster(detectCores() - 1)
  clusterExport(cl, list("simData", "popSim", "f0", "f0se", "stripMissLow", "stripMissMean", "stripMissHigh"))
  clusterEvalQ(cl, library(tidyverse))
  clusterEvalQ(cl, library(extraDistr))
  clusterEvalQ(cl, library(unmarked))
  dataTemp <- parLapply(cl = cl, X = Surveys, fun = function(x) {return(simData(x, popSim, f0, f0se, stripMissLow, stripMissMean, stripMissHigh))})
  stopCluster(cl)
  saveRDS(object = dataTemp, file = paste("output/data/simdata", "_decline", popSimComb[i, "decline"], "_envtlV",
        popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))
  rm(popSim)
  rm(dataTemp)
  gc()

  })

  print(procTime)

}

# fit models
# loop through population combinations
for (i in 1:nrow(popSimComb)) {

  procTime <- system.time({

  # get simulated data
  surveyData <- readRDS(file = paste("output/data/simdata", "_decline", popSimComb[i, "decline"], "_envtlV",
        popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))
  # pick how many replicates - edit out if using all replicates
  surveyData <- lapply(surveyData, FUN = function(x) {x[1:100]})
  expR <-  log((1 - (popSimComb[i, "decline"] / 100)) ^ (1 / genYears))
  cl <- makeCluster(detectCores() - 1)
  clusterExport(cl, list("fitModels", "pops", "expR"))
  clusterEvalQ(cl, library(tidyverse))
  clusterEvalQ(cl, library(unmarked))
  clusterEvalQ(cl, library(binom))
  resultTemp <- parLapply(cl = cl, X = surveyData, fun = function(x) {return(fitModels(x, pops, expR))})
  stopCluster(cl)
  saveRDS(object = resultTemp, file = paste("output/results/result", "_decline", popSimComb[i, "decline"], "_envtlV",
                popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))
  rm(surveyData)
  rm(resultTemp)
  gc()

  })

  print(procTime)

}

# compile results and output rankings

# add row ID to survey combinations
survSimComb <- survSimComb %>% rowid_to_column("ID")

# loop through population simulation combinations
for (i in 1:nrow(popSimComb)) {
  # get type 1 error
  type1 <- readRDS(file = paste("output/results/result", "_decline", 0, "_envtlV",
      popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))
  # get results
  result <- readRDS(file = paste("output/results/result", "_decline", popSimComb[i, "decline"], "_envtlV",
      popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))

  if (popSimComb[i, "decline"] > 0) {

    # compile type 1 errors
    type1Mat <- matrix(unlist(type1), nrow = length(type1), ncol = length(type1[[1]]), byrow = TRUE)
    dimnames(type1Mat) <- list(NULL, names(type1[[1]]))
    type1Mat <- type1Mat %>% as_tibble() %>% rowid_to_column("ID")
    type1Res <- survSimComb %>% left_join(type1Mat, by = "ID") %>% mutate(type1All = region_decline_mean) %>% dplyr::select(ID, type1All)

    # compile results
    resultMat <- matrix(unlist(result), nrow = length(result), ncol = length(result[[1]]), byrow = TRUE)
    dimnames(resultMat) <- list(NULL, names(result[[1]]))
    resultMat <- resultMat %>% as_tibble() %>% rowid_to_column("ID")
    finalRes <- survSimComb %>% left_join(resultMat, by = "ID") %>% mutate(monitBudget = budget, monitYears = monitRep - 1, survInter = monInter, monitStrat = strat, monitMethod = method, type2All = 1 - region_decline_mean,
                      type2GC = 1 - pop1_decline, type2KC = 1 - pop2_decline, type2MB = 1 - pop3_decline, type2NS = 1 - pop4_decline, type2N = 1 - pop5_decline, type2NW = 1 - pop6_decline,
                      type2SW = 1 - pop7_decline) %>%
                      left_join(type1Res, by = "ID") %>%
                      dplyr::select(monitBudget, monitYears, monitStrat, monitMethod, survIntens, survInter, type1All, type2All, type2GC, type2KC, type2MB, type2NS, type2N, type2NW, type2SW) %>%
                      as_tibble() %>% arrange(monitBudget, type2All)

    # save compiled results
    write.csv(x = finalRes, file = paste("output/results/ranking", "_decline", popSimComb[i, "decline"], "_envtlV",
                  popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".csv", sep = ""))
  }
}

# compile monitoring strategy details

# get surveys
Surveys <- apply(survSimComb, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/surveys/survey", "_budget", x["budget"], "_survIntens",
              as.numeric(x["survIntens"]), "_method", x["method"], "_monInter", as.numeric(x["monInter"]), "_strat", x["strat"], "_monitRep",
              as.numeric(x["monitRep"]), ".rds", sep = ""))})

# loop through monitoring strategies
for (i in 1:nrow(survSimComb)) {
  # get output
  Output <- Surveys[[i]] %>% mutate(HABITAT = Quality, LAND_USE = LUCategory, KOALA_DENS = KD_MEAN, NUM_SITES = SITES) %>%
                dplyr::select(FID, POP, HABITAT, LAND_USE, AREA, KOALA_DENS, NUM_SITES, SIZE, TYPE)

  # write output
  write.csv(x = Output, file = paste("output/surveys/monit_strat", "_monitBudget",
                  as.numeric(survSimComb[i, "budget"]), "_survIntens", as.numeric(survSimComb[i, "survIntens"]), "_monitMethod", survSimComb[i, "method"], "_survInter",
                  as.numeric(survSimComb[i, "monInter"]), "_monitStrat", survSimComb[i, "strat"], "_monitYears", as.numeric(survSimComb[i, "monitRep"]) - 1, ".csv", sep = ""))
}

# test adaptive strategy to abandon sites if koalas not seen for a number of surveys

# generate simulated survey data
# get survey information
# here we choose the survey strategies that usually come out as the top two based on Type 2 errors with a budget of 190 (see report) to evaluate the adaptive strategy
survChosen <- survSimComb[c(56, 65),]
Surveys <- apply(survChosen, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/surveys/survey", "_budget", x["budget"], "_survIntens",
              as.numeric(x["survIntens"]), "_method", x["method"], "_monInter", as.numeric(x["monInter"]), "_strat", x["strat"], "_monitRep",
              as.numeric(x["monitRep"]), ".rds", sep = ""))})
# loop through population combinations
for (i in 1:nrow(popSimComb)) {

  # get population simulations
  popSim <- readRDS(file = paste("output/popsims/popsim", "_decline", popSimComb[i, "decline"], "_envtlV",
        popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))

  # generate data for the adaptive strategy
  # loop through adaptation years - here we use between 2 and 5
  for (yearAdapt in 2:5) {

    procTime <- system.time({

    cl <- makeCluster(detectCores() - 1)
    clusterExport(cl, list("simDataAdapt", "popSim", "f0", "f0se", "stripMissLow", "stripMissMean", "stripMissHigh", "yearAdapt"))
    clusterEvalQ(cl, library(tidyverse))
    clusterEvalQ(cl, library(extraDistr))
    clusterEvalQ(cl, library(unmarked))
    clusterEvalQ(cl, library(abind))
    dataTemp <- parLapply(cl = cl, X = Surveys, fun = function(x) {return(simDataAdapt(x, popSim, f0, f0se, stripMissLow, stripMissMean, stripMissHigh, yearAdapt))})
    stopCluster(cl)
    saveRDS(object = dataTemp, file = paste("output/adapt/data/simdata_adapt", "_decline", popSimComb[i, "decline"], "_envtlV",
          popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], "_yearAdapt", yearAdapt, ".rds", sep = ""))
    rm(dataTemp)
    gc()

    })

    print(procTime)

  }

  # generate data for the standard strategy
  procTime <- system.time({

  cl <- makeCluster(detectCores() - 1)
  clusterExport(cl, list("simData", "popSim", "f0", "f0se", "stripMissLow", "stripMissMean", "stripMissHigh"))
  clusterEvalQ(cl, library(tidyverse))
  clusterEvalQ(cl, library(extraDistr))
  clusterEvalQ(cl, library(unmarked))
  dataTemp <- parLapply(cl = cl, X = Surveys, fun = function(x) {return(simData(x, popSim, f0, f0se, stripMissLow, stripMissMean, stripMissHigh))})
  stopCluster(cl)
  saveRDS(object = dataTemp, file = paste("output/adapt/data/simdata_stand", "_decline", popSimComb[i, "decline"], "_envtlV",
        popSimComb[i, "envtlV"], "_envtlCorr", popSimComb[i, "envtlCorr"], ".rds", sep = ""))
  rm(dataTemp)
  gc()

  })

  print(procTime)

  rm(popSim)
  gc()
}

# fit models to adaptive and standard strategies but only for the 0%, 10%, 25% and 50% declines
popChosen <- popSimComb[which((popSimComb[ ,"decline"] == 0) | (popSimComb[ ,"decline"] == 10) | (popSimComb[ ,"decline"] == 25) | (popSimComb[ ,"decline"] == 50)), ]

# loop through population combinations
for (i in 1:nrow(popChosen)) {

  # fit models for the adaptive strategy
  # loop through adaptation years - here we use between 2 and 10
  for (yearAdapt in 2:5) {

    procTime <- system.time({

    # get simulated data for adaptive strategy
    surveyData <- readRDS(file = paste("output/adapt/data/simdata_adapt", "_decline", popChosen[i, "decline"], "_envtlV",
          popChosen[i, "envtlV"], "_envtlCorr", popChosen[i, "envtlCorr"], "_yearAdapt", yearAdapt, ".rds", sep = ""))
    # PICK HOW MANY REPLICATES - edit out if using all replicates
    surveyData <- lapply(surveyData, FUN = function(x) {x[1:100]})
    expR <-  log((1 - (popChosen[i, "decline"] / 100)) ^ (1 / genYears))
    cl <- makeCluster(detectCores() - 1)
    clusterExport(cl, list("fitModels", "pops", "expR"))
    clusterEvalQ(cl, library(tidyverse))
    clusterEvalQ(cl, library(unmarked))
    clusterEvalQ(cl, library(binom))
    resultTemp <- parLapply(cl = cl, X = surveyData, fun = function(x) {return(fitModels(x, pops, expR))})
    stopCluster(cl)
    saveRDS(object = resultTemp, file = paste("output/adapt/results/result_adapt", "_decline", popChosen[i, "decline"], "_envtlV",
                  popChosen[i, "envtlV"], "_envtlCorr", popChosen[i, "envtlCorr"], "_yearAdapt", yearAdapt, ".rds", sep = ""))
    rm(surveyData)
    rm(resultTemp)
    gc()

    })

    print(procTime)
  }

  # fit models for the standard strategy
  procTime <- system.time({

  # get simulated data for adaptive strategy
  surveyData <- readRDS(file = paste("output/adapt/data/simdata_stand", "_decline", popChosen[i, "decline"], "_envtlV",
        popChosen[i, "envtlV"], "_envtlCorr", popChosen[i, "envtlCorr"], ".rds", sep = ""))
  # PICK HOW MANY REPLICATES - edit out if using all replicates
  surveyData <- lapply(surveyData, FUN = function(x) {x[1:100]})
  expR <-  log((1 - (popChosen[i, "decline"] / 100)) ^ (1 / genYears))
  cl <- makeCluster(detectCores() - 1)
  clusterExport(cl, list("fitModels", "pops", "expR"))
  clusterEvalQ(cl, library(tidyverse))
  clusterEvalQ(cl, library(unmarked))
  clusterEvalQ(cl, library(binom))
  resultTemp <- parLapply(cl = cl, X = surveyData, fun = function(x) {return(fitModels(x, pops, expR))})
  stopCluster(cl)
  saveRDS(object = resultTemp, file = paste("output/adapt/results/result_stand", "_decline", popChosen[i, "decline"], "_envtlV",
                popChosen[i, "envtlV"], "_envtlCorr", popChosen[i, "envtlCorr"], ".rds", sep = ""))
  rm(surveyData)
  rm(resultTemp)
  gc()

  })

  print(procTime)
}

# compile adaptive and standard monitoring strategy details for the 1st ranked and 2nd ranked strategy
compiledAdapt_Rank1 <- matrix(NA, nrow = 4, ncol = 6)
dimnames(compiledAdapt_Rank1) <- list(NULL, c("decline", "standard_error", "yearAdapt2_error", "yearAdapt3_error", "yearAdapt4_error",
                        "yearAdapt5_error"))
compiledAdapt_Rank1[ ,"decline"] <- c(0, 10, 20, 50)
compiledAdapt_Rank2 <- matrix(NA, nrow = 4, ncol = 6)
dimnames(compiledAdapt_Rank2) <- list(NULL, c("decline", "standard_error", "yearAdapt2_error", "yearAdapt3_error", "yearAdapt4_error",
                        "yearAdapt5_error"))
compiledAdapt_Rank2[ ,"decline"] <- c(0, 10, 20, 50)

# loop through declines
Declines <- c(0, 10, 25, 50)
for (i in 1:length(Declines)) {

  popChosenTemp <- popChosen[which(popChosen[, "decline"] == Declines[i]), ]

  #get the errors for the standard approach
  meanStand_Rank1 <- mean(apply(popChosenTemp, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/adapt/results/result_stand", "_decline",
              x["decline"], "_envtlV", x["envtlV"], "_envtlCorr", x["envtlCorr"], ".rds", sep = ""))[[1]]["region_decline_mean"]}))
  meanStand_Rank2 <- mean(apply(popChosenTemp, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/adapt/results/result_stand", "_decline",
                          x["decline"], "_envtlV", x["envtlV"], "_envtlCorr", x["envtlCorr"], ".rds", sep = ""))[[2]]["region_decline_mean"]}))
  if (i > 1) {
    meanStand_Rank1 <- 1 - meanStand_Rank1
    meanStand_Rank2 <- 1 - meanStand_Rank2
  }
  compiledAdapt_Rank1[i, "standard_error"] <- meanStand_Rank1
  compiledAdapt_Rank2[i, "standard_error"] <- meanStand_Rank2

  #get the errors for the adaptive approach
  for (j in 2:5) {
    meanAdapt_Rank1 <- mean(apply(popChosenTemp, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/adapt/results/result_adapt", "_decline",
                x["decline"], "_envtlV", x["envtlV"], "_envtlCorr", x["envtlCorr"], "_yearAdapt", j, ".rds", sep = ""))[[1]]["region_decline_mean"]}))
    meanAdapt_Rank2 <- mean(apply(popChosenTemp, MARGIN = 1, FUN = function(x) {readRDS(file = paste("output/adapt/results/result_adapt", "_decline",
                            x["decline"], "_envtlV", x["envtlV"], "_envtlCorr", x["envtlCorr"], "_yearAdapt", j, ".rds", sep = ""))[[2]]["region_decline_mean"]}))
    if (i > 1) {
      meanAdapt_Rank1 <- 1 - meanAdapt_Rank1
      meanAdapt_Rank2 <- 1 - meanAdapt_Rank2
    }
    compiledAdapt_Rank1[i, j + 1] <- meanAdapt_Rank1
    compiledAdapt_Rank2[i, j + 1] <- meanAdapt_Rank2
  }
}

# write outputs
write.csv(x = compiledAdapt_Rank1, file = "output/adapt/results/adapt_compiled_results_rank1.csv")
write.csv(x = compiledAdapt_Rank2, file = "output/adapt/results/adapt_compiled_results_rank2.csv")

max_stars_repo_name: jonrhodes/koalamonitoring

max_issues_repo_name: jonrhodes/koalamonitoring

max_forks_repo_name: jonrhodes/koalamonitoring