tests/testthat/testdata/generate_testdata.R

In danielwilhelm/R-CS-ranks: Statistical Tools for Ranks

library(csranks)
set.seed(100)

####################
### vcov.lmranks ###
####################

for (covariates in c(TRUE,FALSE)) {
  
  # draw data
  n <- 10000
  if (covariates) {
    X <- rnorm(n)
    W <- matrix(rnorm(n*2), n, 2)
    Y <- X + rowSums(W) + rnorm(n,0,0.5)  
    W <- cbind(1,W)  
  } else {
    X <- rnorm(n)
    Y <- X + rnorm(n,0,0.5) 
    W <- matrix(1,n,1)
  }
  
  # compute ranks
  RY <- frank(Y, increasing=TRUE)
  RX <- frank(X, increasing=TRUE)
  
  
  # ------- compute asymptotic variance "by hand"
  
  Ifn <- function(u, v) return( u<=v )
  
  # first stage
  res <- lm(RX ~ W-1)
  Wgammahat <- predict(res)
  nuhat <- resid(res)
  gammahat <- coef(res)
  
  # outcome equation
  res <- lm(RY~RX+W-1)
  rhohat <- coef(res)[1]
  betahat <- coef(res)[-1]
  epsilonhat <- resid(res)
  
  # construct h1
  h1 <- epsilonhat * nuhat
  
  # construct h2
  h2fn <- function(xy) mean((Ifn(xy[2],Y)-rhohat*Ifn(xy[1],X)-c(W%*%betahat)) * nuhat)
  h2 <- apply(cbind(X,Y), 1, h2fn)
  
  # construct h3
  h3fn <- function(x) mean(epsilonhat * (Ifn(x,X)-Wgammahat))
  h3 <- sapply(X, h3fn)
  
  # compute asymptotic variance
  sigma2hat <- mean((h1+h2+h3)^2) / mean(nuhat^2)^2
  
  # save the result
  save(sigma2hat, Y, W, X, n, h1, h2, h3, 
       file = file.path("tests", "testthat", "testdata",
                        paste0("lmranks_cov_sigmahat_covariates_", covariates, ".rda")))
}

for (n in c(10,50,100)) {
  
  # draw data
  
  X <- rnorm(n)
  W <- matrix(rnorm(n*2), n, 2)
  Y <- X + rowSums(W) + rnorm(n,0,0.5)  
  W <- cbind(1,W)  

  
  # compute ranks
  RY <- frank(Y, increasing=TRUE)
  RX <- frank(X, increasing=TRUE)
  
  
  # ------- compute asymptotic variance "by hand"
  
  Ifn <- function(u, v) return( u<=v )
  
  # first stage
  res <- lm(RX ~ W-1)
  Wgammahat <- predict(res)
  nuhat <- resid(res)
  gammahat <- coef(res)
  
  # outcome equation
  res <- lm(RY~RX+W-1)
  rhohat <- coef(res)[1]
  betahat <- coef(res)[-1]
  epsilonhat <- resid(res)
  
  # construct h1
  h1 <- epsilonhat * nuhat
  
  # construct h2
  h2fn <- function(xy) mean((Ifn(xy[2],Y)-rhohat*Ifn(xy[1],X)-c(W%*%betahat)) * nuhat)
  h2 <- apply(cbind(X,Y), 1, h2fn)
  
  # construct h3
  h3fn <- function(x) mean(epsilonhat * (Ifn(x,X)-Wgammahat))
  h3 <- sapply(X, h3fn)
  
  # compute asymptotic variance
  sigma2hat <- mean((h1+h2+h3)^2) / mean(nuhat^2)^2
  
  # save the result
  save(sigma2hat, Y, W, X, n, h1, h2, h3, 
       file = file.path("tests", "testthat", "testdata",
                        paste0("lmranks_cov_sigmahat_n_", n, ".rda")))
}

########################
### increasing=FALSE ###
########################

# draw data
  
X <- rnorm(n)
W <- matrix(rnorm(n*2), n, 2)
Y <- X + rowSums(W) + rnorm(n,0,0.5)  
W <- cbind(1,W)  


# compute ranks
RY <- frank(Y, increasing=FALSE, omega=1)
RX <- frank(X, increasing=FALSE, omega=1)
  
  
# ------- compute asymptotic variance "by hand"
  
Ifn <- function(u, v) return( u>=v )
  
# first stage
res <- lm(RX ~ W-1)
Wgammahat <- predict(res)
nuhat <- resid(res)
gammahat <- coef(res)
  
# outcome equation
res <- lm(RY~RX+W-1)
rhohat <- coef(res)[1]
betahat <- coef(res)[-1]
epsilonhat <- resid(res)

# construct h1
h1 <- epsilonhat * nuhat

# construct h2
h2fn <- function(xy) mean((Ifn(xy[2],Y)-rhohat*Ifn(xy[1],X)-c(W%*%betahat)) * nuhat)
h2 <- apply(cbind(X,Y), 1, h2fn)

# construct h3
h3fn <- function(x) mean(epsilonhat * (Ifn(x,X)-Wgammahat))
h3 <- sapply(X, h3fn)

# compute asymptotic variance
sigma2hat <- mean((h1+h2+h3)^2) / mean(nuhat^2)^2

# save the result
save(sigma2hat, Y, W, X, n, h1, h2, h3, 
     file = file.path("tests", "testthat", "testdata", "lmranks_cov_sigmahat_increasing_FALSE.rda"))

############################
### Hoeffding's variance ###
############################

set.seed(100)

# draw data
n <- 1000
X <- rnorm(n)
Y <- X + rnorm(n,0,0.5) 

# ------- compute Hoeffding's variance "by hand"

n <- length(Y)
rhohat <- cor(frank(Y, omega=1/2, increasing=TRUE),frank(X, omega=1/2, increasing=TRUE))
oX <- outer(X,X,'<=')
oY <- outer(Y,Y,'<=')
FhatX <- colMeans(oX)
FhatY <- colMeans(oY)
FhatAvgX <- function(x) mean(colMeans(outer(X,rep(x,n),'<=')*oY))
FhatAvgY <- function(y) mean(colMeans(outer(Y,rep(y,n),'<=')*oX))
psihatX <- sapply(X, FhatAvgX) - FhatX*mean(FhatY)
psihatY <- sapply(Y, FhatAvgY) - mean(FhatX)*FhatY

W <- 3*( (2*FhatX-1)*(2*FhatY-1) + 4*psihatX + 4*psihatY )

# compute asymptotic variance
sigma2hat <- var(W)

save(sigma2hat, Y, W, X, n, file = file.path("tests", "testthat", "testdata",
                                             "lmranks_cov_sigmahat_Hoefding.rda"))

############################
### vcov.grouped_lmranks ###
############################

set.seed(100)

for (covariates in c(TRUE,FALSE)) {
  # draw data
  n <- 1000
  G <- c(rep(1,n/2),rep(2,n/2))
  beta <- G
  if (covariates) {
    X <- rnorm(n)
    W <- matrix(rnorm(n*2), n, 2)
    Y <- beta/2 + beta*X + rowSums(W) + rnorm(n,0,0.5)  
    W <- cbind(1,W)  
  } else {
    X <- rnorm(n)
    Y <- beta/2 + beta*X + rnorm(n,0,0.5) 
    W <- matrix(1,n,1)
  }
  
  # compute ranks
  RY <- frank(Y, increasing=TRUE, omega=1)
  RX <- frank(X, increasing=TRUE, omega=1)
  RY1 <- RY[1:(n/2)]; RX1 <- RX[1:(n/2)]; W1 <- W[1:(n/2),]
  RY2 <- RY[(n/2+1):n]; RX2 <- RX[(n/2+1):n]; W2 <- W[(n/2+1):n,]
  
  
  # ------- compute asymptotic variance "by hand"
  
  Ifn <- function(u, v) return( u<=v )
  
  # first stage
  res1 <- lm(RX1 ~ W1-1)
  Wgammahat1 <- predict(res1)
  nuhat1 <- resid(res1)
  gammahat1 <- coef(res1)
  res2 <- lm(RX2 ~ W2-1)
  Wgammahat2 <- predict(res2)
  nuhat2 <- resid(res2)
  gammahat2 <- coef(res2)
  nuhat <- c(nuhat1,nuhat2)
  Wgammahat <- c(Wgammahat1,Wgammahat2)
  
  # outcome equation
  res1 <- lm(RY1~RX1+W1-1)
  rhohat1 <- coef(res1)[1]
  betahat1 <- coef(res1)[-1]
  epsilonhat1 <- resid(res1)
  res2 <- lm(RY2~RX2+W2-1)
  rhohat2 <- coef(res2)[1]
  betahat2 <- coef(res2)[-1]
  epsilonhat2 <- resid(res2)
  epsilonhat <- c(epsilonhat1,epsilonhat2)
  
  # construct H1
  H11 <- c(epsilonhat1 * nuhat1, rep(0,n/2))
  H12 <- c(rep(0,n/2), epsilonhat2 * nuhat2)
  
  # construct H2
  H21fn <- function(xy) mean((G==1) * (Ifn(xy[2],Y)-rhohat1*Ifn(xy[1],X)-c(W%*%betahat1)) * nuhat)
  H21 <- apply(cbind(X,Y), 1, H21fn)
  H22fn <- function(xy) mean((G==2) * (Ifn(xy[2],Y)-rhohat2*Ifn(xy[1],X)-c(W%*%betahat2)) * nuhat)
  H22 <- apply(cbind(X,Y), 1, H22fn)
  
  # construct H3
  H31fn <- function(x) mean((G==1) * epsilonhat * (Ifn(x,X)-Wgammahat))
  H31 <- sapply(X, H31fn)
  H32fn <- function(x) mean((G==2) * epsilonhat * (Ifn(x,X)-Wgammahat))
  H32 <- sapply(X, H32fn)
  
  # compute asymptotic variance
  sigma2hat <- c(mean((H11+H21+H31)^2) / mean((G==1)*nuhat^2)^2,
                 mean((H12+H22+H32)^2) / mean((G==2)*nuhat^2)^2)
  
  # save the result
  save(sigma2hat, Y, W, X, G, n, H11, H12, H21, H22, H31, H32,
       file = file.path("tests", "testthat", "testdata",
                                               paste0("grouped_lmranks_cov_sigmahat_covariates_", covariates, ".rda")))
}

for (omega in c(0, 0.4)) {
  # draw data
  n <- 1000
  G <- c(rep(1,n/2),rep(2,n/2))
  beta <- G
  X <- rnorm(n)
  W <- matrix(rnorm(n*2), n, 2)
  Y <- beta/2 + beta*X + rowSums(W) + rnorm(n,0,0.5)  
  W <- cbind(1,W)  
  
  # compute ranks
  RY <- frank(Y, increasing=TRUE, omega=omega)
  RX <- frank(X, increasing=TRUE, omega=omega)
  RY1 <- RY[1:(n/2)]; RX1 <- RX[1:(n/2)]; W1 <- W[1:(n/2),]
  RY2 <- RY[(n/2+1):n]; RX2 <- RX[(n/2+1):n]; W2 <- W[(n/2+1):n,]
  
  
  # ------- compute asymptotic variance "by hand"
  
  Ifn <- function(u, v) return( (u<=v) * omega + (u<v) * (1-omega))
  
  # first stage
  res1 <- lm(RX1 ~ W1-1)
  Wgammahat1 <- predict(res1)
  nuhat1 <- resid(res1)
  gammahat1 <- coef(res1)
  res2 <- lm(RX2 ~ W2-1)
  Wgammahat2 <- predict(res2)
  nuhat2 <- resid(res2)
  gammahat2 <- coef(res2)
  nuhat <- c(nuhat1,nuhat2)
  Wgammahat <- c(Wgammahat1,Wgammahat2)
  
  # outcome equation
  res1 <- lm(RY1~RX1+W1-1)
  rhohat1 <- coef(res1)[1]
  betahat1 <- coef(res1)[-1]
  epsilonhat1 <- resid(res1)
  res2 <- lm(RY2~RX2+W2-1)
  rhohat2 <- coef(res2)[1]
  betahat2 <- coef(res2)[-1]
  epsilonhat2 <- resid(res2)
  epsilonhat <- c(epsilonhat1,epsilonhat2)
  
  # construct H1
  H11 <- c(epsilonhat1 * nuhat1, rep(0,n/2))
  H12 <- c(rep(0,n/2), epsilonhat2 * nuhat2)
  
  # construct H2
  H21fn <- function(xy) mean((G==1) * (Ifn(xy[2],Y)-rhohat1*Ifn(xy[1],X)-c(W%*%betahat1)) * nuhat)
  H21 <- apply(cbind(X,Y), 1, H21fn)
  H22fn <- function(xy) mean((G==2) * (Ifn(xy[2],Y)-rhohat2*Ifn(xy[1],X)-c(W%*%betahat2)) * nuhat)
  H22 <- apply(cbind(X,Y), 1, H22fn)
  
  # construct H3
  H31fn <- function(x) mean((G==1) * epsilonhat * (Ifn(x,X)-Wgammahat))
  H31 <- sapply(X, H31fn)
  H32fn <- function(x) mean((G==2) * epsilonhat * (Ifn(x,X)-Wgammahat))
  H32 <- sapply(X, H32fn)
  
  # compute asymptotic variance
  sigma2hat <- c(mean((H11+H21+H31)^2) / mean((G==1)*nuhat^2)^2,
                 mean((H12+H22+H32)^2) / mean((G==2)*nuhat^2)^2)
  
  # save the result
  save(sigma2hat, Y, W, X, G, n, H11, H12, H21, H22, H31, H32,
       file = file.path("tests", "testthat", "testdata",
                        paste0("grouped_lmranks_cov_sigmahat_omega_", omega, ".rda")))
}

for (n in c(10, 50, 100)) {
  # draw data
  omega <- 0.4
  G <- c(rep(1,n/2),rep(2,n/2))
  beta <- G
  X <- rnorm(n)
  W <- matrix(rnorm(n*2), n, 2)
  Y <- beta/2 + beta*X + rowSums(W) + rnorm(n,0,0.5)  
  W <- cbind(1,W)  
  
  # compute ranks
  RY <- frank(Y, increasing=TRUE, omega=omega)
  RX <- frank(X, increasing=TRUE, omega=omega)
  RY1 <- RY[1:(n/2)]; RX1 <- RX[1:(n/2)]; W1 <- W[1:(n/2),]
  RY2 <- RY[(n/2+1):n]; RX2 <- RX[(n/2+1):n]; W2 <- W[(n/2+1):n,]
  
  
  # ------- compute asymptotic variance "by hand"
  
  Ifn <- function(u, v) return( (u<=v) * omega + (u<v) * (1-omega))
  
  # first stage
  res1 <- lm(RX1 ~ W1-1)
  Wgammahat1 <- predict(res1)
  nuhat1 <- resid(res1)
  gammahat1 <- coef(res1)
  res2 <- lm(RX2 ~ W2-1)
  Wgammahat2 <- predict(res2)
  nuhat2 <- resid(res2)
  gammahat2 <- coef(res2)
  nuhat <- c(nuhat1,nuhat2)
  Wgammahat <- c(Wgammahat1,Wgammahat2)
  
  # outcome equation
  res1 <- lm(RY1~RX1+W1-1)
  rhohat1 <- coef(res1)[1]
  betahat1 <- coef(res1)[-1]
  epsilonhat1 <- resid(res1)
  res2 <- lm(RY2~RX2+W2-1)
  rhohat2 <- coef(res2)[1]
  betahat2 <- coef(res2)[-1]
  epsilonhat2 <- resid(res2)
  epsilonhat <- c(epsilonhat1,epsilonhat2)
  
  # construct H1
  H11 <- c(epsilonhat1 * nuhat1, rep(0,n/2))
  H12 <- c(rep(0,n/2), epsilonhat2 * nuhat2)
  
  # construct H2
  H21fn <- function(xy) mean((G==1) * (Ifn(xy[2],Y)-rhohat1*Ifn(xy[1],X)-c(W%*%betahat1)) * nuhat)
  H21 <- apply(cbind(X,Y), 1, H21fn)
  H22fn <- function(xy) mean((G==2) * (Ifn(xy[2],Y)-rhohat2*Ifn(xy[1],X)-c(W%*%betahat2)) * nuhat)
  H22 <- apply(cbind(X,Y), 1, H22fn)
  
  # construct H3
  H31fn <- function(x) mean((G==1) * epsilonhat * (Ifn(x,X)-Wgammahat))
  H31 <- sapply(X, H31fn)
  H32fn <- function(x) mean((G==2) * epsilonhat * (Ifn(x,X)-Wgammahat))
  H32 <- sapply(X, H32fn)
  
  # compute asymptotic variance
  sigma2hat <- c(mean((H11+H21+H31)^2) / mean((G==1)*nuhat^2)^2,
                 mean((H12+H22+H32)^2) / mean((G==2)*nuhat^2)^2)
  
  # save the result
  save(sigma2hat, Y, W, X, G, n, H11, H12, H21, H22, H31, H32,
       file = file.path("tests", "testthat", "testdata",
                        paste0("grouped_lmranks_cov_sigmahat_n_", n, ".rda")))
}