ich knabbere seit geraumer Zeit an einem Problem, für das ich noch keine Antwort finden konnte. Und zwar will ich den Code unten ausführen. Dieser simuliert mir 1000 mal einen Illness-Death Prozess https://journals.plos.org/plosone/artic ... ne.0123489 und ich möchte daraus eine Übergangswahrscheinlichkeit auf zwei Arten schätzen (Details sind erst mal nicht so wichtig). Dazu berechne ich eine empirische Übergangsmatrix mitdem Befehl etm aus dem etm-package. Jedoch stürzt RStudio mir jedes mal (nach einer gewissen Anzahl Iterationen, so 10-50) ab mit folgender Fehlermeldung :
error:: Mat::operator() index out of bounds
terminate called after throwing an instance of 'std::logic_error'
what(): Mat::operator() index out of bounds
This application has requested the Runtime to terminate it in an unusual way.
Please contact the application's support team for more information.
RStudio zeigt mir die Fehlermeldung '<p>R encountered a fatal error.</p>The session was terminated.' Danach muss ich die Sitzung neu starten usw.
Code: Alles auswählen
rm(list = ls())
library(etm)
library(mstate)
#markovian case
BSiter <- 300
bs_one_iter <- function(n, data, max_id, times=x, h, tra){
multinom <- sample(1:n, replace=T)
data_bs <- subset(data, id==data$id[multinom[1]])
for(k in 2:n){
data_temp <- subset(data, id==data$id[multinom[k]])
data_temp$id <- data_temp$id + max_id * k
data_bs <- rbind(data_bs, data_temp)
}
etmAJE <- etm(data_bs, state.names = 0:2, tra=tra, cens.name=3, s=s, t=t, covariance=FALSE, delta.na=FALSE)
PijAJE <- c(0,etmAJE$est[i + 1,j + 1,])
timesAJE <- c(s,etmAJE$time)
good_ind <- which(data$entry <= s & data$exit >s & data$from==i)
rel_ids <- unique(data$id[good_ind])
npeople <- length(rel_ids)
data_s <- data[data$id %in% rel_ids,]
#Remove all intervals ending before s
data_s <- data_s[data_s$exit >= s,]
#Truncate all intervals starting before s
data_s$entry[data_s$entry <s]<-s
etm_s <- etm(data_s, state.names = 0:2, tra=tra, cens.name=3, s=s, t=t, covariance=FALSE, delta.na=FALSE)
return(k)
}
for (k in 1:1000){
# simulation
n <- 100
id <- 1:n
rate_0 <- 1/15
start_state <- rep(0,n)
T_0 <- rexp(n,rate_0) # time in state 0
p <- 0.7
r <- runif(n,0,1)
state_1 <- 1*(r < p) + 2*(p <= r & r < p+(1-p)/2) + 3*(1-(1-p)/2 <= r)
rate_1 <- 1/30
T_1 <- rexp(n,rate_1) # time in state 1
q <- 0.8
state_2 <- rbinom(n,1,1-q) + 2
# results in time-continuous markovian multi-state process'
# state 3 = censored
# next we transform the generated data into a data matrix that shows transition times and states per patient
all_datas <- matrix(rep(0,5*2*n),nrow = 2*n)
count <- 1
for (i in id){
if(state_1[i] == 3 | state_1[i] == 2){
all_datas[count,] <- c(i,0,state_1[i],0,T_0[i])
}
else{
all_datas[count,] <- c(i,0,1,0,T_0[i])
count <- count + 1
all_datas[count,] <- c(i,1,state_2[i],T_0[i],T_0[i] + T_1[i])
}
count <- count + 1
}
datas <- data.frame(all_datas[-which(all_datas[,1] == 0),])
names(datas) <- c("id","from","to","entry","exit")
tra <- matrix(c(F,T,T,F,F,T,F,F,F),3,3,byrow=T)
# from illness to death
i <- 1
j <- 2
# time stamp
s <- 5
t <- 100
# estimate transition probabilities
etmAJE_markov <- etm(datas, state.names = 0:2, tra=tra, cens.name=3, s=s, t=t, covariance=FALSE, delta.na=FALSE)
PijAJE_markov <- c(0,etmAJE_markov$est[i + 1,j + 1,])
timesAJE_markov <- c(s,etmAJE_markov$time)
good_ind_markov <- which(datas$entry <= s & datas$exit >s & datas$from == i)
rel_ids_markov <- unique(datas$id[good_ind_markov])
npeople_markov <- length(rel_ids_markov)
# "good" data
datas_s <- data.frame(datas[datas$id %in% rel_ids_markov,])
# Remove all intervals ending before s
datas_s <- datas_s[datas_s$exit >= s,]
#Truncate all intervals starting before s
datas_s$entry[datas_s$entry < s] <- s
etm_s_markov <- etm(datas_s, state.names = 0:2, tra=tra, cens.name=3, s=s, t=t, covariance=FALSE, delta.na=FALSE)
F0_markov <- c(0,etm_s_markov$est[i + 1,j + 1,])
timesTE_markov <- c(s,etm_s_markov$time)
# Computing the Aalen Johansen and Titman estimators
h <- 0.5
x_markov <- seq(s,t,h)
# Computing the test statistic
max_id_markov <- max(datas$id)
BSstat_pred <- replicate(BSiter, bs_one_iter(n, datas, max_id_markov, times=x_markov, h, tra))
}
Hat jemand eine Idee, wo genau der Fehler passieren könnte und was an den simulaierten Daten das Problem ist?
Btw. Wenn ich die Daten wie folgt simuliere, läuft der Prozess problemlos durch.
Ich habe auch schon versucht den Autor des packages zu kontaktieren, jedoch stimmt seine Adresse in der Dokumentation nicht mehr... Die Konsole stürzt mir bei der Ausführung ebenso ab.
Code: Alles auswählen
# zuerst wie oben bis Kommentar ' # simulation'
n <- 1000
id <- 1:n
size <- 50
box <- rexp(size,1/2)
lambda_0 <- 10
nc0 <- rpois(n,lambda_0) #number of coins to draw in state 0
r <- runif(n,0,1)
p <- 0.7
state_1 <- 1*(r < p) + 2*(p <= r & r < p+(1-p)/2) + 3*(1-(1-p)/2 <= r)
lambda_1 <- 15
nc1 <- rpois(n,lambda_1) #number of coins to draw in state 1
q <- 0.8
state_2 <- rbinom(n,1,1-q) + 2
T_0 <- rep(0,n)
T_1 <- rep(0,n)
# 'draw' coins from the box without replacement and note the sum to get the time of transfer between states
# results in a time-continuous non-markov multi state process
for (i in id){
zero_ind <- sample(size,nc0[i])
zero_samp <- box[zero_ind]
T_0[i] <- sum(zero_samp)
remain <- box[-zero_samp]
one_ind <- sample(size - nc0[i],nc1[i])
one_samp <- remain[one_ind]
T_1[i] <- sum(one_samp)
}
... # dann weiter wie ab dem Kommentar '# next we transform the generated data into a data matrix that shows transition times and states per patient'
Entschuldigt, dass der Code relativ lang ist, aber leider habe ich ziemlich viele Abhängigkeiten darin, die nicht so ohne weiteres auszumerzen sind (es handelt sich schon um ein vereinfachtes Beispiel. Deswegen gibt die Funktion bs_one_iter auch Unsinn zurück.). Bei Fragen helfe ich natürlich.
Danke und viele Grüße
Tonelock
Du müsstest herausfinden, bei welcher Datenkonstellation der Absturz stattfindet 
