CODEWITHSUNDEEP
rfor-loopmatchinglogistic-regressionsurvival-analysis
lukeg
lukeg

Reputation: 1357

writing a program to introduce permutations

Basically I want to write a program that will randomise the order of my data n times, then complete a survival analysis and plot the output over n

So lets take the following generic data from the matching() package and create a dataset of treated and non-treated people. Link to package

set.seed(123)

library(Matching)
data(lalonde)

lalonde$age_cat <- with(lalonde, ifelse(age < 24, 1, 2))
attach(lalonde)

lalonde$ID <- 1:length(lalonde$age)


#The covariates we want to match on
X = cbind(age_cat, educ, black, hisp, married, nodegr, u74, u75, re75, re74)
#The covariates we want to obtain balance on
BalanceMat <- cbind(age_cat, educ, black, hisp, married, nodegr, u74, u75, re75, re74,
                    I(re74*re75))
genout <- GenMatch(Tr=treat, X=X, BalanceMatrix=BalanceMat, estimand="ATE", M=1,
                   pop.size=16, max.generations=10, wait.generations=1)
detach(lalonde)

# now lets pair the the non-treated collisions to the treated
# BUT lets pair WITHOUT REPLACEMENT

mout <- Match(Y=NULL, Tr=lalonde$treat, X=X,
              Weight.matrix=genout, M=2,
              replace=FALSE, ties=TRUE)

summary(mout)
# we see that for 130 treated observations, we have 260 non-treated
# this is because we set M=2
# and yes length(lalonde$age[lalonde$treat==0]) == 260 but just follow me please
# but this was done for a specific reason

# now lets create a table for our 130+260 collisions
treated <- lalonde[mout$index.treated,]
# now we only want one occurence of the treated variables
library(dplyr)
treat_clean <- treated %>%
  group_by(ID) %>%
  slice(1)

non.treated <- lalonde[mout$index.control,]

# finally we can combine to form one clear data.set
matched.data <- rbind(treat_clean, non.treated)

We can now do a conditional logistic regression to determine the OR associated with re78 (money earned in 1987) and treatment. For this we need the survival package. Link to package

library(survival)

Lets say a success occurs if the occupant earns more than 8125 in 1978

matched.data$success <- with(matched.data, ifelse(re78 > 8125, 1, 0))

output <- clogit(success ~ treat, matched.data, method = 'efron')

summary(output)

so we see that the OR for the treated (treat=1) is 1.495

We can save this as:

iteration.1 <- exp(output$coefficients[1])

Now we read from the matching package (link) that for replace = FALSE Note that if FALSE, the order of matches generally matters. Matches will be found in the same order as the data are sorted

So what I want to do create a function that will for n times

Is essenece I want to introduce permutations into the analysis. How can this be done, when lets say n=5

Upvotes: 0

Views: 117

Answers (2)

Rorschach
Rorschach

Reputation: 32446

You can use sample to introduce permutations

data(lalonde)
lalonde$age_cat <- with(lalonde, ifelse(age < 24, 1, 2))
lalonde$ID <- 1:length(lalonde$age)
n <- 5
res <- rep(NA, n)
for (i in 1:n) {
    lalonde <- lalonde[sample(1:nrow(lalonde)), ] # randomise order
    ## rest of code 
    res[i] <- exp(output$coefficients[1])
}

plot(1:n, res, main="Odds Ratios")

Upvotes: 1

r2evans
r2evans

Reputation: 160607

I'm a big fan of replicate for something like this:

X <- cbind(...)                         # what you had before
BalanceMat <- cbind(...)                # ditto
lalonde$ID <- seq.int(nrow(lalonde))

results <- replicate(1000, {
    ## not certain if it's just $ID order that matters
    lalonde$ID <- sample(nrow(lalonde))
    ## lalonde <- lalonde[ sample(nrow(lalonde)), ]

    ## ...
    ## rest of your computation
    ## ...

    #### optionally return everything 
    ## output
    #### return just the minimum
    exp(output$coefficients[1])
})

#### if you returned output earlier, you'll need this, otherwise not
## coef <- exp(sapply(results, function(z) z$coefficients[1]))

## plot as needed

I don't know if you meant just the order of ID matters or if the order of the entire database; adjust the first couple lines of the replicate loop accordingly.

Upvotes: 1

Related Questions