Кривые мощности с использованием оценок мощности из симстуды - PullRequest
Новая
       69

Кривые мощности с использованием оценок мощности из симстуды

0 голосов
/ 04 мая 2020

На слайдах 8-10 следующего представления: http://www.columbia.edu/~cjd11/charles_dimaggio/DIRE/resources/R/power.pdf, мощность моделируется для t-теста с двумя образцами. Затем строится кривая мощности с использованием кода на слайде 14 с использованием ggplot: 

# Power curve for sample size
sizes <- seq(100, 3000, by=50)
powerVals <- vector(length=length(sizes))
for (i in 1:length(sizes)) {
powerVals[i] <- fnPower(intercept=1,beta=0.15,sd=1,nTotalSubjects=sizes[i],nTreated=sizes[i]/2,
nIterations=100)$power
}
powerCalc <- data.frame(sizes, powerVals)
ggplot(powerCalc) + aes(x=sizes, y=powerVals) + geom_smooth()

Можно ли построить кривые мощности, используя оценки мощности, которые были получены из пакета simstudy? Как бы я go сделал это? Вот мой код для оценки мощности с использованием simstudy для справки: 

library(simstudy)
library(nlme)
library(gendata)
library(data.table)
library(geepack)

set.seed(12345)

clusterDef <- defDataAdd(varname = "u_3", dist = "normal", formula = 0, variance = 25.77) #cluster-level random effect
patError <- defDataAdd(varname = "error", dist = "normal", formula = 0, variance = 38.35) #error term 

#Generate cluster-level data
cohortsw <- genData(3, id = "cluster")
cohortsw <- addColumns(clusterDef, cohortsw) 
cohortswTm <- addPeriods(cohortsw, nPeriods = 6, idvars = "cluster", perName = "period")
cohortstep <- trtStepWedge(cohortswTm, "cluster", nWaves = 3, lenWaves = 1, startPer = 1, grpName = "Ijt")
cohortstep

#Generate individual patient-level data
pat <- genCluster(cohortswTm, cLevelVar = "timeID", numIndsVar = 5, level1ID = "id")
pat

dx <- merge(pat[, .(cluster, period, id)], cohortstep, by = c("cluster", "period")) 
dx <- addColumns(patError, dx)

setkey(dx, id, cluster, period)

#Define outcome y 
outDef <- defDataAdd(varname = "y", formula = "17.87 + 5.0*Ijt - 5.42*I(period == 1) - 5.72*I(period == 2) - 7.03*I(period == 3) - 6.13*I(period == 4) - 9.13*I(period == 5) + u_3 + error", dist = "normal")
dx <- addColumns(outDef, dx)

#Fit GLMM model to simulated dataset
model1 <- lme(y ~ factor(period) + factor(Ijt), random = ~1|cluster, data = dx, method = "REML")
summary(model1)

#Power analysis 
powercrosssw <- function(nclus = 3, clsize = 5) {

  cohortsw <- genData(nclus, id = "cluster")
  cohortsw <- addColumns(clusterDef, cohortsw) 
  cohortswTm <- addPeriods(cohortsw, nPeriods = 6, idvars = "cluster", perName = "period")
  cohortstep <- trtStepWedge(cohortswTm, "cluster", nWaves = 3, lenWaves = 1, startPer = 1, grpName = "Ijt")

  pat <- genCluster(cohortswTm, cLevelVar = "timeID", numIndsVar = clsize, level1ID = "id")

  dx <- merge(pat[, .(cluster, period, id)], cohortstep, by = c("cluster", "period"))
  dx <- addColumns(patError, dx)

  setkey(dx, id, cluster, period)

  return(dx)

}

bresult <- NULL
presult <- NULL
eresult <- NULL

intercept <- NULL 
trt <- NULL 
timecoeff1 <- NULL
timecoeff2 <- NULL
timecoeff3 <- NULL
timecoeff4 <- NULL
timecoeff5 <- NULL
ranclus <- NULL
error <- NULL 

i=1

while (i < 1000) {

  cohortsw <- powercrosssw()

  #Fit multi-level model to simulated dataset
  model1 <- tryCatch(lme(y ~ factor(period) + factor(Ijt), data = dx, random = ~1|cluster, method = "REML"),
                     warning = function(w) { "warning" }
  )

  if (! is.character(model1)) {

    coeff <- coef(summary(model1))["factor(Ijt)1", "Value"]
    pvalue <- coef(summary(model1))["factor(Ijt)1", "p-value"]
    error <- coef(summary(model1))["factor(Ijt)1", "Std.Error"]
    bresult <- c(bresult, coeff)
    presult <- c(presult, pvalue)
    eresult <- c(eresult, error)

    i <- i + 1
  }
}
...