Nothing
tests/testthat/test-7-koch-et-al.R
In getspanel: General-to-Specific Modelling of Panel Data
test_that("Equal Breaks as in Koch et al are identified",{
skip_on_cran()
# Prepare Data
set.seed(1230)
#data <- read.csv("CO2DriversEU_dataset.csv")
data("EU_emissions_road")
data <- EU_emissions_road
data$lgdp_sq <- data$lgdp^2
data$transport.emissions_pc <- data$transport.emissions/data$pop
data$ltransport.emissions_pc <- log(data$transport.emissions_pc)
data$L1.ltransport.emissions_pc <- NA
# For each country, shift the values of 'ltransport.emissions_pc' by one position
for (i in unique(data$country)) {
# Extract the 'ltransport.emissions_pc' values for the current country
current_country_values <- data$ltransport.emissions_pc[data$country == i]
# Shift the values by one position and insert an NA value at the beginning
shifted_values <- c(NA, current_country_values[-length(current_country_values)])
# Assign the shifted values to the corresponding rows in 'L1.ltransport.emissions_pc'
data$L1.ltransport.emissions_pc[data$country == i] <- shifted_values
}
# data %>%
# dplyr::group_by(country) %>%
# dplyr::mutate(test = dplyr::lag(ltransport.emissions_pc)) %>%
# dplyr::ungroup() %>%
# dplyr::mutate(istrue = test == L1.ltransport.emissions_pc) %>%
# dplyr::distinct(istrue)
# Group specification
EU15 <- c("Austria", "Belgium", "Germany", "Denmark", "Spain", "Finland",
"France", "United Kingdom", "Ireland", "Italy", "Luxembourg",
"Netherlands", "Greece", "Portugal", "Sweden")
# EU16 <- c("Croatia", "Bulgaria", "Cyprus","Czech Republic", "Estonia",
# "Hungary", "Lithuania", "Latvia", "Malta", "Poland", "Romania",
# "Slovak Republic", "Slovenia", "Switzerland", "Iceland",
# "Norway")
# EU31 <- c(EU15, EU16)
###### Analysis:
# Prepare sample and data
sample <- EU15
dat <- data[data$country %in% sample & data$year>=1995,]
# Break analysis:
is1 <- isatpanel(
data = dat,
formula = ltransport.emissions_pc ~ lgdp + lgdp_sq + lpop,
index = c("country", "year"),
effect = "twoways",
iis = T,
fesis = T,
t.pval=.05,
print.searchinfo = FALSE
)
# Break analysis:
is2 <- isatpanel(
data = dat,
formula = ltransport.emissions_pc ~ lgdp + lgdp_sq + lpop,
index = c("country", "year"),
effect = "twoways",
iis = T,
fesis = T,
t.pval=.01,
print.searchinfo = FALSE
)
# Break analysis:
is3 <- isatpanel(
data = dat,
formula = ltransport.emissions_pc ~ lgdp + lgdp_sq + lpop,
index = c("country", "year"),
effect = "twoways",
iis = T,
fesis = T,
t.pval=.001,
print.searchinfo = FALSE
)
get_countries <- function(x){
df <- x$estimateddata
indicators <- x$isatpanel.result$aux$mX
indicators <- indicators[,!colnames(indicators) %in% names(df)]
df <- cbind(df,indicators)
identify_indicator_timings(df, uis_breaks = NULL)
}
is1_coef <- coef(is1$isatpanel.result)[get_countries(is1)$fesis$name]
is2_coef <- coef(is2$isatpanel.result)[get_countries(is2)$fesis$name]
is3_coef <- coef(is3$isatpanel.result)[get_countries(is3)$fesis$name]
is1_coef <- is1_coef[order(names(is1_coef))]
is2_coef <- is2_coef[order(names(is2_coef))]
is3_coef <- is3_coef[order(names(is3_coef))]
is1_tib <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisIreland.2015",
"fesisLuxembourg.2007",
"fesisSweden.2001"),
coef = c(-0.103,
-0.105,
-0.087,
-0.148,
-0.136,
-0.095))
is2_tib <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2015",
"fesisSweden.2001"),
coef = c(-0.123,
-0.131,
-0.192,
-0.103))
is3_tib <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisLuxembourg.2015",
"fesisSweden.2001"),
coef = c(-0.128,
-0.108,
-0.127,
-0.214,
-0.110))
# Check that Equal Breaks as in Koch et al are identified
is1_tib_estimated <- data.frame(name = names(is1_coef[is1_coef < 0]),
coef = as.numeric(round(is1_coef[is1_coef < 0], 3)))
expect_equal(is1_tib, is1_tib_estimated)
is2_tib_estimated <- data.frame(name = names(is2_coef[is2_coef < 0]),
coef = as.numeric(round(is2_coef[is2_coef < 0], 3)))
expect_equal(is2_tib, is2_tib_estimated)
is3_tib_estimated <- data.frame(name = names(is3_coef[is3_coef < 0]),
coef = as.numeric(round(is3_coef[is3_coef < 0], 3)))
expect_equal(is3_tib, is3_tib_estimated)
# Check that Equal Break Uncertainties as in Koch et al are identified"
break_is1 <- break_uncertainty( is1, interval = 0.99)
break_is2 <- break_uncertainty( is2, interval = 0.99)
break_is3 <- break_uncertainty( is3, interval = 0.99)
break_is1 <- break_is1[break_is1$coef < 0,]
break_is2 <- break_is2[break_is2$coef < 0,]
break_is3 <- break_is3[break_is3$coef < 0,]
row.names(break_is1) <- 1:nrow(break_is1)
row.names(break_is2) <- 1:nrow(break_is2)
row.names(break_is3) <- 1:nrow(break_is3)
is1_tib_break <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisIreland.2015",
"fesisLuxembourg.2007",
"fesisSweden.2001"),
tci = c(2,
2,
3,
1,
1,
2))
is2_tib_break <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2015",
"fesisSweden.2001"),
tci = c(2,
1,
1,
2))
is3_tib_break <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisLuxembourg.2015",
"fesisSweden.2001"),
tci = c(2,
3,
2,
1,
3))
expect_equal(is1_tib_break, break_is1[,c("name", "tci")])
expect_equal(is2_tib_break, break_is2[,c("name", "tci")])
expect_equal(is3_tib_break, break_is3[,c("name", "tci")])
# Check that Koch et al: Standard Error Corrections are correct
is1_robust_nocluster <- robust_isatpanel(is1, cluster = FALSE)$robust
is2_robust_nocluster <- robust_isatpanel(is2, cluster = FALSE)$robust
is3_robust_nocluster <- robust_isatpanel(is3, cluster = FALSE)$robust
is1_robust_nocluster_round <- round(is1_robust_nocluster[is1_robust_nocluster[,1,drop = FALSE]<0 & grepl("^fesis",row.names(is1_robust_nocluster)),2,drop = FALSE],3)
is2_robust_nocluster_round <- round(is2_robust_nocluster[is2_robust_nocluster[,1,drop = FALSE]<0 & grepl("^fesis",row.names(is2_robust_nocluster)),2,drop = FALSE],3)
is3_robust_nocluster_round <- round(is3_robust_nocluster[is3_robust_nocluster[,1,drop = FALSE]<0 & grepl("^fesis",row.names(is3_robust_nocluster)),2,drop = FALSE],3)
is1_robust_nocluster_ready <- data.frame(name = row.names(is1_robust_nocluster_round), se = is1_robust_nocluster_round[,1])
is2_robust_nocluster_ready <- data.frame(name = row.names(is2_robust_nocluster_round), se = is2_robust_nocluster_round[,1])
is3_robust_nocluster_ready <- data.frame(name = row.names(is3_robust_nocluster_round), se = is3_robust_nocluster_round[,1])
is1_robust_nocluster_ready <- is1_robust_nocluster_ready[order(is1_robust_nocluster_ready$name),]
is2_robust_nocluster_ready <- is2_robust_nocluster_ready[order(is2_robust_nocluster_ready$name),]
is3_robust_nocluster_ready <- is3_robust_nocluster_ready[order(is3_robust_nocluster_ready$name),]
row.names(is1_robust_nocluster_ready) <- NULL
row.names(is2_robust_nocluster_ready) <- NULL
row.names(is3_robust_nocluster_ready) <- NULL
is1_robust_nocluster_paper <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisIreland.2015",
"fesisLuxembourg.2007",
"fesisSweden.2001"),
se = c(0.006,
0.008,
0.012,
0.031,
0.012,
0.008))
is2_robust_nocluster_paper <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2015",
"fesisSweden.2001"),
se = c(0.014,
0.011,
0.032,
0.01))
is3_robust_nocluster_paper <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisLuxembourg.2015",
"fesisSweden.2001"),
se = c(0.018,
0.02,
0.016,
0.034,
0.015))
expect_equal(is1_robust_nocluster_paper, is1_robust_nocluster_ready)
expect_equal(is2_robust_nocluster_paper, is2_robust_nocluster_ready)
expect_equal(is3_robust_nocluster_paper, is3_robust_nocluster_ready)
})
# ##### uncertainty for break dates
#
# ####This file compiles the uncertainty functions needed for sis
#
#
# library(mvtnorm)
#
# ### Function for rotation
# rotate <- function(x) t(apply(x, 2, rev))
#
# ### Function to extract matrices
# mextrc <- function(x, p=0, k=0){
# # k=2*(plus minus), i.e. total range, and p is where it starts relative to the center, so k is the size, and p is the starting point
#
# mid <- NROW(x)/2
# out <- 2*k
# rotate <- function(x) t(apply(x, 2, rev))
# subx <- x[(mid+p-k):(mid+p-k+out+1),(mid+p-k):(mid+p-k+out+1)]
# subx <- rotate(rotate(rotate(x)))[(mid+p-k):(mid+p-k+out+1),(mid+p-k):(mid+p-k+out+1)]
# output <- rotate(rotate(rotate(subx)))
# return(output)
# }
#
# #####################################################
#
# ####### Function to compute approximate SIS uncertainty
#
# isattime <- function(m, delta, plot=FALSE){
#
# #######################################
#
# lim <- 2*m
# sigma <- diag(2*m)
# a_lim <- m-1
#
# ########-----------Constructing general matrix
#
# up0 <- (matrix(NA, nrow=2*m, ncol=1))
# up <- (matrix(NA, nrow=2*m, ncol=(m)))
# upg.a <- seq(1:lim)
# upg <- c(rev(upg.a*(-1)), upg.a)
# upg[(lim+1):(lim+m-1)]
#
#
# ###loop to generate
# for (j in 1:m){
#
# #T0 values
# up0[j, 1] <- delta/2*((m+1)-j)
# up0[j+m, 1] <- delta/2*j
#
# #Tp+m values
# ##loop over rows
# for (l in 1:(m+1))
# {
# #j is the column
# up[l,j] <- delta/2 *(((m+1)-l)-j)
# up[l+m-1, j] <- delta/2 * upg[(lim-1-j+l)]
# }
# }
#
#
# up.tot <- cbind(up0, up)
#
# #######------Constructing the general A Matrix
#
# A <- matrix(0, nrow=lim, ncol=lim)
# for (i in 1:lim){
# A[i, 1:i] <- 1
# }
# A
# Ar <- rotate(A)
# Arm <- (-1)*rotate(rotate(rotate(A)))
# Az <- matrix(0, nrow=lim, ncol=lim)
# Ag <- rbind(cbind(Az, Ar), cbind(Arm, Az))
# Ag
# Agb <- Ag*(-1)
#
# ###A0
# A0 <- mextrc(Agb, 0, a_lim )*(-1)
# A0gsigma <- A0 %*% sigma %*% t(A0)
# p0 <- pmvnorm(mean=rep(0, NROW(A0gsigma)), sigma=A0gsigma, lower=rep(-Inf, NROW(A0gsigma)), upper=up.tot[,1] )[1]
# p0
#
# p_m <- matrix(NA, m, 1)
#
# ### up to m
# for (i in 1:m){
# Ai <- mextrc(Agb, i, a_lim )
# Aisigma <- Ai %*% sigma %*% t(Ai)
# pi <- pmvnorm(mean=rep(0, NROW(Aisigma)), sigma=Aisigma, lower=rep(-Inf, NROW(Aisigma)), upper=up.tot[,(i+1)] )[1]
# p_m[i,1] <- pi
# }
#
# p.tot <- as.matrix(c(rev(p_m), p0, p_m))
# p.tot
# colSums(p.tot)
#
# x <- seq(from = -m, to=m, by=1)
#
# if (plot){
# plot(x, p.tot, type='hist')
# }
# return(p.tot)
# }
#
# #######################################################
# #################################################
#
#
#
# ############# Input
#
# model.names <- c("EU15_levelpc_p0.05", "EU15_levelpc_p0.01", "EU15_levelpc_p0.001")
#
# dat.input <- list() #break magnitudes (estimated coefficients)
# breaknames <- list() #break title for overview
# se <- c() #s.e. of regression (i.e. sqrt of estimated error variance)
#
# # Extract coeffs and se from model .Rds:
# for(k in 1:length(model.names)){
# #mod <- model.names[k]
#
# #raw <- readRDS(paste0("../../A Break Detection/Analysis4a (level)/Analysis4a/", mod, ".Rds"))
# raw <- list(is1, is2, is3)[k][[1]]
#
# finalmodel <- raw$isatpanel.result$mean.results
#
# coeffs <- finalmodel$coef[stringr::str_detect(row.names(finalmodel), "fesis")]
# breaks <- row.names(finalmodel)[stringr::str_detect(row.names(finalmodel), "fesis")] %>%
# stringr::str_remove("fesis")
#
# sigm <- sqrt(raw$isatpanel.result$sigma2)
#
#
# dat.input[[k]] <- coeffs
# se[k] <- sigm
# breaknames[[k]] <- breaks
# }
#
#
# ############ Calibration parameters
# interval <- 0.99 #approx. level of interval. CI level will be at least > interval. So 0.99 is a 99% CI, and interval will be the integer that gets at least > 99% coverage
# m <- 15 #maximum range of interval (can leave at 15 for default).
#
# set.seed(1234)
#
# #######################################
# #################
#
# for(k in 1:length(model.names)){
# dat <- dat.input[[k]]/se[k] #standardizing break magnitudes
# prob <- matrix(NA, NROW(dat), 2*m+1)
# ci <- matrix(NA, NROW(dat), 4)
#
# for (i in 1:NROW(dat)){
# #i <- 1
# time <- isattime(m, abs(dat[i])) #external function to compute approx. CI
#
# prob[i,] <- t(time)
# cu.p <- matrix(NA, (length(time)+1)/2, 1)
# for (j in 1:((length(time)+1)/2))
# {
# start <- (length(time)+1)/2
# cu.p[j,1] <- sum(time[((start-j+1):(start+j-1))])
# }
#
# thresh <- interval
# ci.v <- min(which(abs(cu.p)>thresh))
# ci.d <- ci.v - 1
# cu.p[ci.v]
#
# ci[i,1] <- ci.d #interval is given by estimated date + - ci.d
# ci[i,2] <- cu.p[ci.v] #coverage probability
# ci[i,3] <- breaknames[[k]][i] #which break is considered
# ci[i,4] <- dat.input[[k]][i] #break magnitude
#
# }
#
# colnames(ci) <- c("half_range", "coverage", "break_title", "break_magnitude")
#
# cat(paste0("\n Model ", k, "/", length(model.names), " done."))
# }
Try the getspanel package in your browser
Any scripts or data that you put into this service are public.
getspanel documentation built on June 8, 2025, 12:51 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
test_that("Equal Breaks as in Koch et al are identified",{
skip_on_cran()
# Prepare Data
set.seed(1230)
#data <- read.csv("CO2DriversEU_dataset.csv")
data("EU_emissions_road")
data <- EU_emissions_road
data$lgdp_sq <- data$lgdp^2
data$transport.emissions_pc <- data$transport.emissions/data$pop
data$ltransport.emissions_pc <- log(data$transport.emissions_pc)
data$L1.ltransport.emissions_pc <- NA
# For each country, shift the values of 'ltransport.emissions_pc' by one position
for (i in unique(data$country)) {
# Extract the 'ltransport.emissions_pc' values for the current country
current_country_values <- data$ltransport.emissions_pc[data$country == i]
# Shift the values by one position and insert an NA value at the beginning
shifted_values <- c(NA, current_country_values[-length(current_country_values)])
# Assign the shifted values to the corresponding rows in 'L1.ltransport.emissions_pc'
data$L1.ltransport.emissions_pc[data$country == i] <- shifted_values
}
# data %>%
# dplyr::group_by(country) %>%
# dplyr::mutate(test = dplyr::lag(ltransport.emissions_pc)) %>%
# dplyr::ungroup() %>%
# dplyr::mutate(istrue = test == L1.ltransport.emissions_pc) %>%
# dplyr::distinct(istrue)
# Group specification
EU15 <- c("Austria", "Belgium", "Germany", "Denmark", "Spain", "Finland",
"France", "United Kingdom", "Ireland", "Italy", "Luxembourg",
"Netherlands", "Greece", "Portugal", "Sweden")
# EU16 <- c("Croatia", "Bulgaria", "Cyprus","Czech Republic", "Estonia",
# "Hungary", "Lithuania", "Latvia", "Malta", "Poland", "Romania",
# "Slovak Republic", "Slovenia", "Switzerland", "Iceland",
# "Norway")
# EU31 <- c(EU15, EU16)
###### Analysis:
# Prepare sample and data
sample <- EU15
dat <- data[data$country %in% sample & data$year>=1995,]
# Break analysis:
is1 <- isatpanel(
data = dat,
formula = ltransport.emissions_pc ~ lgdp + lgdp_sq + lpop,
index = c("country", "year"),
effect = "twoways",
iis = T,
fesis = T,
t.pval=.05,
print.searchinfo = FALSE
)
# Break analysis:
is2 <- isatpanel(
data = dat,
formula = ltransport.emissions_pc ~ lgdp + lgdp_sq + lpop,
index = c("country", "year"),
effect = "twoways",
iis = T,
fesis = T,
t.pval=.01,
print.searchinfo = FALSE
)
# Break analysis:
is3 <- isatpanel(
data = dat,
formula = ltransport.emissions_pc ~ lgdp + lgdp_sq + lpop,
index = c("country", "year"),
effect = "twoways",
iis = T,
fesis = T,
t.pval=.001,
print.searchinfo = FALSE
)
get_countries <- function(x){
df <- x$estimateddata
indicators <- x$isatpanel.result$aux$mX
indicators <- indicators[,!colnames(indicators) %in% names(df)]
df <- cbind(df,indicators)
identify_indicator_timings(df, uis_breaks = NULL)
}
is1_coef <- coef(is1$isatpanel.result)[get_countries(is1)$fesis$name]
is2_coef <- coef(is2$isatpanel.result)[get_countries(is2)$fesis$name]
is3_coef <- coef(is3$isatpanel.result)[get_countries(is3)$fesis$name]
is1_coef <- is1_coef[order(names(is1_coef))]
is2_coef <- is2_coef[order(names(is2_coef))]
is3_coef <- is3_coef[order(names(is3_coef))]
is1_tib <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisIreland.2015",
"fesisLuxembourg.2007",
"fesisSweden.2001"),
coef = c(-0.103,
-0.105,
-0.087,
-0.148,
-0.136,
-0.095))
is2_tib <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2015",
"fesisSweden.2001"),
coef = c(-0.123,
-0.131,
-0.192,
-0.103))
is3_tib <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisLuxembourg.2015",
"fesisSweden.2001"),
coef = c(-0.128,
-0.108,
-0.127,
-0.214,
-0.110))
# Check that Equal Breaks as in Koch et al are identified
is1_tib_estimated <- data.frame(name = names(is1_coef[is1_coef < 0]),
coef = as.numeric(round(is1_coef[is1_coef < 0], 3)))
expect_equal(is1_tib, is1_tib_estimated)
is2_tib_estimated <- data.frame(name = names(is2_coef[is2_coef < 0]),
coef = as.numeric(round(is2_coef[is2_coef < 0], 3)))
expect_equal(is2_tib, is2_tib_estimated)
is3_tib_estimated <- data.frame(name = names(is3_coef[is3_coef < 0]),
coef = as.numeric(round(is3_coef[is3_coef < 0], 3)))
expect_equal(is3_tib, is3_tib_estimated)
# Check that Equal Break Uncertainties as in Koch et al are identified"
break_is1 <- break_uncertainty( is1, interval = 0.99)
break_is2 <- break_uncertainty( is2, interval = 0.99)
break_is3 <- break_uncertainty( is3, interval = 0.99)
break_is1 <- break_is1[break_is1$coef < 0,]
break_is2 <- break_is2[break_is2$coef < 0,]
break_is3 <- break_is3[break_is3$coef < 0,]
row.names(break_is1) <- 1:nrow(break_is1)
row.names(break_is2) <- 1:nrow(break_is2)
row.names(break_is3) <- 1:nrow(break_is3)
is1_tib_break <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisIreland.2015",
"fesisLuxembourg.2007",
"fesisSweden.2001"),
tci = c(2,
2,
3,
1,
1,
2))
is2_tib_break <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2015",
"fesisSweden.2001"),
tci = c(2,
1,
1,
2))
is3_tib_break <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisLuxembourg.2015",
"fesisSweden.2001"),
tci = c(2,
3,
2,
1,
3))
expect_equal(is1_tib_break, break_is1[,c("name", "tci")])
expect_equal(is2_tib_break, break_is2[,c("name", "tci")])
expect_equal(is3_tib_break, break_is3[,c("name", "tci")])
# Check that Koch et al: Standard Error Corrections are correct
is1_robust_nocluster <- robust_isatpanel(is1, cluster = FALSE)$robust
is2_robust_nocluster <- robust_isatpanel(is2, cluster = FALSE)$robust
is3_robust_nocluster <- robust_isatpanel(is3, cluster = FALSE)$robust
is1_robust_nocluster_round <- round(is1_robust_nocluster[is1_robust_nocluster[,1,drop = FALSE]<0 & grepl("^fesis",row.names(is1_robust_nocluster)),2,drop = FALSE],3)
is2_robust_nocluster_round <- round(is2_robust_nocluster[is2_robust_nocluster[,1,drop = FALSE]<0 & grepl("^fesis",row.names(is2_robust_nocluster)),2,drop = FALSE],3)
is3_robust_nocluster_round <- round(is3_robust_nocluster[is3_robust_nocluster[,1,drop = FALSE]<0 & grepl("^fesis",row.names(is3_robust_nocluster)),2,drop = FALSE],3)
is1_robust_nocluster_ready <- data.frame(name = row.names(is1_robust_nocluster_round), se = is1_robust_nocluster_round[,1])
is2_robust_nocluster_ready <- data.frame(name = row.names(is2_robust_nocluster_round), se = is2_robust_nocluster_round[,1])
is3_robust_nocluster_ready <- data.frame(name = row.names(is3_robust_nocluster_round), se = is3_robust_nocluster_round[,1])
is1_robust_nocluster_ready <- is1_robust_nocluster_ready[order(is1_robust_nocluster_ready$name),]
is2_robust_nocluster_ready <- is2_robust_nocluster_ready[order(is2_robust_nocluster_ready$name),]
is3_robust_nocluster_ready <- is3_robust_nocluster_ready[order(is3_robust_nocluster_ready$name),]
row.names(is1_robust_nocluster_ready) <- NULL
row.names(is2_robust_nocluster_ready) <- NULL
row.names(is3_robust_nocluster_ready) <- NULL
is1_robust_nocluster_paper <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisIreland.2015",
"fesisLuxembourg.2007",
"fesisSweden.2001"),
se = c(0.006,
0.008,
0.012,
0.031,
0.012,
0.008))
is2_robust_nocluster_paper <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2015",
"fesisSweden.2001"),
se = c(0.014,
0.011,
0.032,
0.01))
is3_robust_nocluster_paper <- data.frame(name = c("fesisFinland.2000",
"fesisGermany.2002",
"fesisIreland.2011",
"fesisLuxembourg.2015",
"fesisSweden.2001"),
se = c(0.018,
0.02,
0.016,
0.034,
0.015))
expect_equal(is1_robust_nocluster_paper, is1_robust_nocluster_ready)
expect_equal(is2_robust_nocluster_paper, is2_robust_nocluster_ready)
expect_equal(is3_robust_nocluster_paper, is3_robust_nocluster_ready)
})
# ##### uncertainty for break dates
#
# ####This file compiles the uncertainty functions needed for sis
#
#
# library(mvtnorm)
#
# ### Function for rotation
# rotate <- function(x) t(apply(x, 2, rev))
#
# ### Function to extract matrices
# mextrc <- function(x, p=0, k=0){
# # k=2*(plus minus), i.e. total range, and p is where it starts relative to the center, so k is the size, and p is the starting point
#
# mid <- NROW(x)/2
# out <- 2*k
# rotate <- function(x) t(apply(x, 2, rev))
# subx <- x[(mid+p-k):(mid+p-k+out+1),(mid+p-k):(mid+p-k+out+1)]
# subx <- rotate(rotate(rotate(x)))[(mid+p-k):(mid+p-k+out+1),(mid+p-k):(mid+p-k+out+1)]
# output <- rotate(rotate(rotate(subx)))
# return(output)
# }
#
# #####################################################
#
# ####### Function to compute approximate SIS uncertainty
#
# isattime <- function(m, delta, plot=FALSE){
#
# #######################################
#
# lim <- 2*m
# sigma <- diag(2*m)
# a_lim <- m-1
#
# ########-----------Constructing general matrix
#
# up0 <- (matrix(NA, nrow=2*m, ncol=1))
# up <- (matrix(NA, nrow=2*m, ncol=(m)))
# upg.a <- seq(1:lim)
# upg <- c(rev(upg.a*(-1)), upg.a)
# upg[(lim+1):(lim+m-1)]
#
#
# ###loop to generate
# for (j in 1:m){
#
# #T0 values
# up0[j, 1] <- delta/2*((m+1)-j)
# up0[j+m, 1] <- delta/2*j
#
# #Tp+m values
# ##loop over rows
# for (l in 1:(m+1))
# {
# #j is the column
# up[l,j] <- delta/2 *(((m+1)-l)-j)
# up[l+m-1, j] <- delta/2 * upg[(lim-1-j+l)]
# }
# }
#
#
# up.tot <- cbind(up0, up)
#
# #######------Constructing the general A Matrix
#
# A <- matrix(0, nrow=lim, ncol=lim)
# for (i in 1:lim){
# A[i, 1:i] <- 1
# }
# A
# Ar <- rotate(A)
# Arm <- (-1)*rotate(rotate(rotate(A)))
# Az <- matrix(0, nrow=lim, ncol=lim)
# Ag <- rbind(cbind(Az, Ar), cbind(Arm, Az))
# Ag
# Agb <- Ag*(-1)
#
# ###A0
# A0 <- mextrc(Agb, 0, a_lim )*(-1)
# A0gsigma <- A0 %*% sigma %*% t(A0)
# p0 <- pmvnorm(mean=rep(0, NROW(A0gsigma)), sigma=A0gsigma, lower=rep(-Inf, NROW(A0gsigma)), upper=up.tot[,1] )[1]
# p0
#
# p_m <- matrix(NA, m, 1)
#
# ### up to m
# for (i in 1:m){
# Ai <- mextrc(Agb, i, a_lim )
# Aisigma <- Ai %*% sigma %*% t(Ai)
# pi <- pmvnorm(mean=rep(0, NROW(Aisigma)), sigma=Aisigma, lower=rep(-Inf, NROW(Aisigma)), upper=up.tot[,(i+1)] )[1]
# p_m[i,1] <- pi
# }
#
# p.tot <- as.matrix(c(rev(p_m), p0, p_m))
# p.tot
# colSums(p.tot)
#
# x <- seq(from = -m, to=m, by=1)
#
# if (plot){
# plot(x, p.tot, type='hist')
# }
# return(p.tot)
# }
#
# #######################################################
# #################################################
#
#
#
# ############# Input
#
# model.names <- c("EU15_levelpc_p0.05", "EU15_levelpc_p0.01", "EU15_levelpc_p0.001")
#
# dat.input <- list() #break magnitudes (estimated coefficients)
# breaknames <- list() #break title for overview
# se <- c() #s.e. of regression (i.e. sqrt of estimated error variance)
#
# # Extract coeffs and se from model .Rds:
# for(k in 1:length(model.names)){
# #mod <- model.names[k]
#
# #raw <- readRDS(paste0("../../A Break Detection/Analysis4a (level)/Analysis4a/", mod, ".Rds"))
# raw <- list(is1, is2, is3)[k][[1]]
#
# finalmodel <- raw$isatpanel.result$mean.results
#
# coeffs <- finalmodel$coef[stringr::str_detect(row.names(finalmodel), "fesis")]
# breaks <- row.names(finalmodel)[stringr::str_detect(row.names(finalmodel), "fesis")] %>%
# stringr::str_remove("fesis")
#
# sigm <- sqrt(raw$isatpanel.result$sigma2)
#
#
# dat.input[[k]] <- coeffs
# se[k] <- sigm
# breaknames[[k]] <- breaks
# }
#
#
# ############ Calibration parameters
# interval <- 0.99 #approx. level of interval. CI level will be at least > interval. So 0.99 is a 99% CI, and interval will be the integer that gets at least > 99% coverage
# m <- 15 #maximum range of interval (can leave at 15 for default).
#
# set.seed(1234)
#
# #######################################
# #################
#
# for(k in 1:length(model.names)){
# dat <- dat.input[[k]]/se[k] #standardizing break magnitudes
# prob <- matrix(NA, NROW(dat), 2*m+1)
# ci <- matrix(NA, NROW(dat), 4)
#
# for (i in 1:NROW(dat)){
# #i <- 1
# time <- isattime(m, abs(dat[i])) #external function to compute approx. CI
#
# prob[i,] <- t(time)
# cu.p <- matrix(NA, (length(time)+1)/2, 1)
# for (j in 1:((length(time)+1)/2))
# {
# start <- (length(time)+1)/2
# cu.p[j,1] <- sum(time[((start-j+1):(start+j-1))])
# }
#
# thresh <- interval
# ci.v <- min(which(abs(cu.p)>thresh))
# ci.d <- ci.v - 1
# cu.p[ci.v]
#
# ci[i,1] <- ci.d #interval is given by estimated date + - ci.d
# ci[i,2] <- cu.p[ci.v] #coverage probability
# ci[i,3] <- breaknames[[k]][i] #which break is considered
# ci[i,4] <- dat.input[[k]][i] #break magnitude
#
# }
#
# colnames(ci) <- c("half_range", "coverage", "break_title", "break_magnitude")
#
# cat(paste0("\n Model ", k, "/", length(model.names), " done."))
# }
Try the getspanel package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.