Commit 9e9272cf by Jan Wijffels

Change comparison graphs to functionality from the citizenair package

parent 705d8345
......@@ -7,6 +7,7 @@
#' Defaults to 'citizen'
#' @return an object of class \code{citizenair_comparison_ircel} which is a list with elements
#' \enumerate{
#' \item{valid: logical indicating if we have enough data to do the regression}
#' \item{data: the data containing columns time, REF, SENSOR, OLSQ, lower_CI, upper_CI, lower_PI, upper_PI and e}
#' \item{model: an object of class lm represendint the model SENSOR ~ REF}
#' \item{overview: a data.frame with columns statistic, time, REF, SENSOR, OLSQ, lower_CI, upper_CI, lower_PI, upper_PI and e with the
......@@ -20,6 +21,18 @@
#' value = rnorm(365*2))
#' comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
#' plot(comparison)
#'
#' x <- data.frame(
#' time = rep(seq.POSIXt(from = Sys.time()-364, to = Sys.time(), by = "sec"), 2),
#' group = sample(c("official", "citizen"), size = 365*2, replace = TRUE),
#' value = rnorm(365*2))
#' comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
#'
#' x <- data.frame(
#' time = as.Date(character()),
#' group = character(),
#' value = numeric())
#' comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
citizenair_compare <- function(x, reference = "official", sensor = "citizen"){
# r cmd check happiness
.SD <- OLSQ <- REF <- SENSOR <- lower_CI <- lower_PI <- upper_CI <- upper_PI <- NULL
......@@ -32,32 +45,8 @@ citizenair_compare <- function(x, reference = "official", sensor = "citizen"){
x$group <- txt_recode(x$group, from = c(reference, sensor), to = c("REF", "SENSOR"))
x <- data.table::setDT(x)
mydata <- data.table::dcast.data.table(data = x, formula = time ~ group,
fun.aggregate = mean, na.rm=TRUE,
value.var = "value")
mydata <- data.table::setDF(mydata)
linear_model <- lm(SENSOR ~ REF, data = mydata)
PI <- predict(linear_model, interval = "predict", level = 0.95, newdata = mydata, na.action = na.pass)
PI <- as.data.frame(PI)
colnames(PI) <- c("OLSQ", "lower PI", "upper PI")
CI <- predict(linear_model, interval = "confidence", level = 0.95, newdata = mydata, na.action = na.pass)
CI <- as.data.frame(CI)
colnames(CI) <- c("OLSQ", "lower CI", "upper CI")
mydata[["OLSQ"]] <- CI[["OLSQ"]]
mydata[["lower_CI"]] <- CI[["lower CI"]]
mydata[["upper_CI"]] <- CI[["upper CI"]]
mydata[["lower_PI"]] <- PI[["lower PI"]]
mydata[["upper_PI"]] <- PI[["upper PI"]]
mydata$OLSQ <- CI$OLSQ
mydata$OLSQ <- CI$OLSQ
mydata$OLSQ <- CI$OLSQ
mydata$OLSQ <- CI$OLSQ
mydata$OLSQ <- CI$OLSQ
mydata$e <- mydata$OLSQ - mydata$SENSOR
mydata <- mydata[order(mydata$time, decreasing = FALSE), ]
overview <- data.table::as.data.table(mydata)
overview_statistics <- function(x){
overview <- data.table::as.data.table(x)
overview <- overview[, lapply(.SD, FUN=function(x){
summarystats <- summary(x)
if(inherits(x, c("Date", "POSIXct", "POSIXlt"))){
......@@ -76,8 +65,47 @@ citizenair_compare <- function(x, reference = "official", sensor = "citizen"){
overview <- data.table::setDF(overview)
overview[["statistic"]] <- c("min", "25th Perc", "median", "mean", "75th Perc", "max", "\u0023 Valid", "\u0023 NA's")
overview <- overview[, c("statistic", setdiff(colnames(overview), "statistic"))]
overview
}
if(nrow(x) == 0){
mydata <- data.frame(time = as.Date(character()), SENSOR = numeric(), REF = numeric(),
OLSQ = numeric(),
lower_CI = numeric(),
upper_CI = numeric(),
lower_PI = numeric(),
upper_PI = numeric(),
stringsAsFactors = FALSE)
result <- list(valid = FALSE,
data = mydata,
overview = overview_statistics(mydata),
model = NULL)
class(result) <- "citizenair_comparison_ircel"
return(result)
}else{
mydata <- data.table::dcast.data.table(data = x, formula = time ~ group,
fun.aggregate = mean, na.rm=TRUE,
value.var = "value")
mydata <- data.table::setDF(mydata)
}
linear_model <- lm(SENSOR ~ REF, data = mydata)
PI <- predict(linear_model, interval = "predict", level = 0.95, newdata = mydata, na.action = na.pass)
PI <- as.data.frame(PI)
colnames(PI) <- c("OLSQ", "lower PI", "upper PI")
CI <- predict(linear_model, interval = "confidence", level = 0.95, newdata = mydata, na.action = na.pass)
CI <- as.data.frame(CI)
colnames(CI) <- c("OLSQ", "lower CI", "upper CI")
mydata[["OLSQ"]] <- CI[["OLSQ"]]
mydata[["lower_CI"]] <- CI[["lower CI"]]
mydata[["upper_CI"]] <- CI[["upper CI"]]
mydata[["lower_PI"]] <- PI[["lower PI"]]
mydata[["upper_PI"]] <- PI[["upper PI"]]
mydata$e <- mydata$OLSQ - mydata$SENSOR
mydata <- mydata[order(mydata$time, decreasing = FALSE), ]
overview <- overview_statistics(mydata)
result <- list(data = mydata,
result <- list(valid = TRUE,
data = mydata,
overview = overview,
model = linear_model)
class(result) <- "citizenair_comparison_ircel"
......@@ -97,10 +125,20 @@ citizenair_compare <- function(x, reference = "official", sensor = "citizen"){
#' value = rnorm(365*2))
#' comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
#' plot(comparison)
#'
#' x <- data.frame(
#' time = as.Date(character()),
#' group = character(),
#' value = numeric())
#' comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
#' plot(comparison)
plot.citizenair_comparison_ircel <- function(x, ...){
OLSQ <- REF <- SENSOR <- lower_CI <- lower_PI <- upper_CI <- upper_PI <- NULL
linear_model <- x$model
mydata <- x$data
if(!x$valid){
plot_PI_CI <- ggplot2::ggplot(data = mydata)
}else{
plot_PI_CI <- ggplot2::ggplot(data = mydata, ggplot2::aes(x = REF, y = SENSOR)) +
ggplot2::geom_ribbon(ggplot2::aes(ymin = lower_CI, ymax=upper_CI, fill = "CI"), alpha = 0.6) +
ggplot2::geom_ribbon(ggplot2::aes(ymin = lower_PI, ymax=upper_PI, fill = "PI"), alpha = 0.4) +
......@@ -117,6 +155,7 @@ plot.citizenair_comparison_ircel <- function(x, ...){
"PI"= "grey")) +
ggplot2::scale_colour_manual(name="", values = c("y=x"="red",
"OLSQ"="black"))
}
plot_PI_CI
}
......
......@@ -161,6 +161,7 @@ reactive({
input$uiInput_aggregationlevel
input$uiInput_sheet
input$uiInput_period
input$uiInput_dayofweek
#showNotification("Herberekening wordt uitgevoerd")
#input$uiInput_phenomena
......@@ -172,69 +173,54 @@ reactive({
x <- subset(x, type %in% "official" | (type %in% "citizen" & timeseries_id %in% input$uiInput_sheet))
x <- subset(x, phenomena_id %in% input$uiInput_phenomena)
x <- x[as.Date(x$time) >= (as.Date(input$uiInput_period[1])) & as.Date(x$time) <= (as.Date(input$uiInput_period[2])), ]
if(input$uiInput_dayofweek %in% c("weekdays", "weekend")){
if(input$uiInput_dayofweek %in% "weekdays"){
x <- x[format(as.Date(x$time), "%u") %in% c(1, 2, 3, 4, 5), ]
}else if(input$uiInput_dayofweek %in% "weekend"){
x <- x[format(as.Date(x$time), "%u") %in% c(6, 7), ]
}
}
x <- APPDATA$getCombinedTimeseries(x, type = input$uiInput_aggregationlevel, limit = TRUE)
#output$uiOutput_timeseries <- renderDygraph(plot_timeseries(x))
try({
mydata <- subset(x, phenomena_id %in% head(input$uiInput_phenomena, 1))
if(all(c("citizen", "official") %in% mydata$type) & input$uiInput_aggregationlevel %in% c("min", "hour", "day", "week", "month")){
## DATE/REF/SENSOR
mydata <- x
mydata$group <- ifelse(mydata$type %in% "official", "REF", "SENSOR")
mydata$date <- mydata$time
mydata <- dcast.data.table(data = mydata, formula = date ~ group, fun.aggregate = mean, na.rm=TRUE, value.var = "value")
#mydata <- data.frame(date = as.POSIXct(seq.Date(Sys.Date(), Sys.Date()-999, by = "-1 day"), tz = "UTC"),
# REF = rnorm(1000),
# SENSOR = rnorm(1000))
linear_model <- lm(mydata$SENSOR ~ mydata$REF)
PI <- data.frame(predict(linear_model, interval="predict", level = 0.95))
colnames(PI) <- c("OLSQ", "lower PI", "upper PI")
CI <- data.frame(predict(linear_model, interval="confidence", level = 0.95))
colnames(CI) <- c("OLSQ", "lower CI", "upper CI")
mydata <- cbind(mydata, "OLSQ" = CI$OLSQ,
"lower_CI" = CI$`lower CI`,
"upper_CI"=CI$`upper CI`,
"lower_PI" = PI$`lower PI`,
"upper_PI" = PI$`upper PI`)
mydata$e <- mydata$OLSQ - mydata$SENSOR
mydata <- mydata[order(mydata[,1]), ]
summary <- data.frame(do.call(cbind, lapply(mydata, summary)))[-7,]
summary$date <- as.character(as.POSIXct(summary$date, origin="1970-01-01"))
summary <- rbind(summary,"# Valid"=colSums(!is.na(mydata)))
summary <- rbind(summary,"# NA's"=colSums(is.na(mydata)))
variables <- c("min", "25th Perc", "median", "mean","75th Perc", "max","# Valid", "# NA's")
summary <- cbind("statistiek"=variables,summary)
# outputs of dataframe & summaries
comparison_phenomena <- head(input$uiInput_phenomena, 1)
mydata <- subset(x, phenomena_id %in% comparison_phenomena)
comparison_phenomena <- APPDATA$getPhenomenaLabel(comparison_phenomena)
if(input$uiInput_aggregationlevel %in% c("min", "hour", "day", "week", "month") && nrow(mydata) > 0 && all(c("citizen", "official") %in% mydata$type)){
mydata <- data.frame(time = mydata$time,
group = ifelse(mydata$type %in% "citizen", "SENSOR", ifelse(mydata$type %in% "official", "REFERENCE", NA)),
value = mydata$value,
stringsAsFactors = FALSE)
comparison <- try(citizenair_compare(x = mydata, reference = "REFERENCE", sensor = "SENSOR"))
if(inherits(comparison, "try-error")){
showNotification(sprintf("Er loopt iets fout bij de vergelijking %s", error_message(comparison)),
type = "error", duration = 15)
output$data <- renderTable(NULL)
output$summaryRaw <- renderTable(NULL)
output$uiOutput_comparison <- renderPlot(NULL)
output$uiOutput_comparison_measurement <- renderUI(tags$h4("We analyseren hier uw meting van ", comparison_phenomena, " ... probleem met de regressie"))
}else{
output$uiOutput_comparison_measurement <- renderUI(tags$h4("We analyseren hier uw meting van ", comparison_phenomena))
output$data <- renderTable({
mydata <- timeAverage(mydata, avg.time=input$uiInput_aggregationlevel, data.thresh=85, statistic="mean")
mydata$date <- as.character(as.POSIXct(mydata$date, origin="1970-01-01"))
comparison$data <- data.table::setnames(comparison$data, old = "time", new = "date")
mydata <- timeAverage(comparison$data,
avg.time=input$uiInput_aggregationlevel, data.thresh = 85, statistic = "mean")
mydata$date <- as.character(as.POSIXct(mydata$date, origin = "1970-01-01"))
mydata
})
output$summaryRaw <- renderTable(summary)
# Scatterplot of linear regression + CI & PI
output$uiOutput_comparison <- renderPlot({
mydata <- mydata
linear_model <- linear_model
plot_PI_CI <- ggplot(data = mydata, aes(x = REF, y = SENSOR)) +
geom_ribbon(aes(ymin = lower_CI, ymax=upper_CI, fill = "CI"),alpha = 0.6) +
geom_ribbon(aes(ymin = lower_PI, ymax=upper_PI, fill = "PI"), alpha = 0.4) +
geom_point() +
geom_line(aes(y = OLSQ,col = "OLSQ")) +
geom_line(aes(y = REF,col = "y=x")) +
labs(title = paste("Scatterplot with regression estimates: intercept =",
format(summary(linear_model)$coefficients[1,1], digits = 3),
",slope =",format(summary(linear_model)$coefficients[2,1], digits = 3),
",R² =",format(summary(linear_model)$r.squared, digits = 3)),
x = "Reference value",
y = "Sensor") +
scale_fill_manual(name="",values=c("CI"="blue","PI"= "grey")) +
scale_colour_manual(name="", values = c("y=x"="red","OLSQ"="black"))
#plots[["scatter"]] <- plot_PI_CI
plot_PI_CI
})
output$summaryRaw <- renderTable(comparison$overview)
output$uiOutput_comparison <- renderPlot(plot(comparison))
}
}else{
if(length(comparison_phenomena) > 0){
output$uiOutput_comparison_measurement <- renderUI(tags$h4("We analyseren hier uw meting van ", comparison_phenomena, " - We hebben geen meting van zowel u als een officieel station. Selecteer een andere meetwaarde die u wilt vergelijken en zorg ervoor dat de vergelijkingsbasis per uur/dag of maand is."))
}else{
output$uiOutput_comparison_measurement <- renderUI(tags$h4("We hebben geen meting van zowel u als een officieel station. Selecteer een andere meetwaarde die u wilt vergelijken, haal data op van een officieel meetstation en zorg ervoor dat de vergelijkingsbasis per uur/dag of maand is."))
}
output$data <- renderTable(NULL)
output$summaryRaw <- renderTable(NULL)
output$uiOutput_comparison <- renderPlot(NULL)
}
})
})
```
......@@ -260,7 +246,8 @@ fillCol(
tags$div(
tags$h4("De Vlaamse Milieu Maatschappij heeft meetresultaten die aan de strengste kwaliteitseisen voldoen. Hieronder tonen
we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van jouw. Verander ook de VMM-meetplaats eens, een andere meetplaats komt misschien beter met de jouwe overeen!"),
tags$h5("Onderstaande grafiek toont alle meetwaardes van uw toestel uitgezet tov alle meetwaardes van een toestel van de Vlaamse Milieu Maatschappij")),
tags$h5("Onderstaande grafiek toont alle meetwaardes van uw toestel uitgezet tov alle meetwaardes van een toestel van de Vlaamse Milieu Maatschappij. Uw toestel moet in lijn liggen met een officieel toestel. Dit kunt u zien door te kijken of uw meetpunten op 1 lijn liggen en heel dicht liggen tegen de rode lijn.")),
htmlOutput(outputId = "uiOutput_comparison_measurement"),
fillRow(dropdownButton(
tags$h3("Selectie dagen van de week"),
radioButtons(inputId = "uiInput_dayofweek", label = "Vergelijk op basis van",
......@@ -273,7 +260,7 @@ we een eenvoudige statistische vergelijking tussen onze meetgegevens en die van
plotOutput("uiOutput_comparison", height = "100%"),
flex = c(NA, 2)
),
flex = c(NA, 2)
flex = c(NA, NA, 2)
)
```
......
......@@ -18,6 +18,7 @@ Defaults to 'citizen'}
\value{
an object of class \code{citizenair_comparison_ircel} which is a list with elements
\enumerate{
\item{valid: logical indicating if we have enough data to do the regression}
\item{data: the data containing columns time, REF, SENSOR, OLSQ, lower_CI, upper_CI, lower_PI, upper_PI and e}
\item{model: an object of class lm represendint the model SENSOR ~ REF}
\item{overview: a data.frame with columns statistic, time, REF, SENSOR, OLSQ, lower_CI, upper_CI, lower_PI, upper_PI and e with the
......@@ -34,4 +35,16 @@ x <- data.frame(
value = rnorm(365*2))
comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
plot(comparison)
x <- data.frame(
time = rep(seq.POSIXt(from = Sys.time()-364, to = Sys.time(), by = "sec"), 2),
group = sample(c("official", "citizen"), size = 365*2, replace = TRUE),
value = rnorm(365*2))
comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
x <- data.frame(
time = as.Date(character()),
group = character(),
value = numeric())
comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
}
......@@ -22,4 +22,11 @@ x <- data.frame(
value = rnorm(365*2))
comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
plot(comparison)
x <- data.frame(
time = as.Date(character()),
group = character(),
value = numeric())
comparison <- citizenair_compare(x, reference = "official", sensor = "citizen")
plot(comparison)
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment