Discussion :

[andy18]

Bonjour,

J'espère vraiment que quelqu'un va pouvoir m'aider, je travaille sur des simulations de sous-échantillonage pour des données de survie et voici mon problème : j'utilise une fonction "sim" qui appelle une fonction "arms" (crée 2 groupes de survie) dans une boucle. Lorsque je lance la fonction "arms" hors de la fonction, elle marche très bien. Mais lorsque je lance la fonction la variable à laquelle j'attribue les résultst de la fonction "arms" n'est pas crée. Voici le code de la fonction "arms" :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
 
 
library(splines)
library(survival)
 
# r = percentage of patients with treatment 1 (between 0 and 1)
# N = number of patients
 
arms <- function(r,N){
 
	y1 <- 0.9 * rexp(ceiling(r*N), rate = 1/4)
	y2 <- 1.1 * rexp(ceiling((1-r)*N), rate = 1/2)
 
	time	<- c(y1, y2)
	tmax <- max(time)
	n1 <- length(y1)
	n2 <- length(y2)
 
	treat1 <- rep(0, n1)
	treat2 <- rep(1, n2)
	treat <- c(treat1, treat2)
 
	censtimes1 <- 16 * runif(n1)			# censoring time vector
	ztimes1 <- pmin( y1, censtimes1 )			# choose which comes first
	status1 <- as.numeric ( censtimes1 > y1 )		# related indicator
 
	censtimes2 <- 11 * runif(n2)			# censoring time vector
	ztimes2 <- pmin( y2, censtimes2 )			# choose which comes first
	status2 <- as.numeric ( censtimes2 > y2 )		# related indicator
 
	status <- c(status1, status2)
 
	samp <- cbind(time, status, treat)
	samp <- data.frame(samp)
 
	fit <- survfit(Surv(time, status) ~ treat, data = samp)
	fitnocens <- survfit(Surv(time, status) ~ treat, data = samp[samp$status==1,])
	#plot(fit, col = c(1,3))
	#legend("topright",legend = c("T1","T2"), col = c("black","green"), pch = 15)
 
	l <- list(y1=y1,y2=y2,tmax=tmax,n1=n1,n2=n2, fit=fit, time=time, status=status, treat=treat, samp=samp, fitnocens=fitnocens)
	return(l)
}

et la fonction "sim", l'étape qui me bloque est en rouge (au début de la fonction)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

#nb = number of simulations
#N = number of patients in the sample
#r = percentage (between 0 and 1) of patients in the 1st treatment
#p1 and p2 = percentages (between 0 and 1) of patients in each treatment in the subsample out of the sample
#subsamp = choose the subsamp1 : "subsample" or 2 : "halfsubsample"
#cens = 0 if no censored data are required in the subsample, 1 otherwise, 
#this option only works for the second subsample (but is required in the "sim" function)
sim <- function(nb, N, r, p1, p2, subsamp, cens){

	coxsample <- rep(NA, nb)
	coxsubsample <- rep(NA, nb)

	for (i in 1:nb){

		mainsample <- arms(r,N)
		subsample <- subsamp( p1, p2, cens)	

		coxsample[i] <- summary(coxph(Surv(mainsample$time, mainsample$status) ~ mainsample$treat))[[7]][[1]] 
		coxsubsample[i] <- summary(coxph(Surv(subsample$time, subsample$status) ~ subsample$treat))[[7]][[1]]
	}

	diff <- coxsample - coxsubsample
	coeffs <- cbind(coxsample ,coxsubsample,diff)

	densample <- density(coeffs[,1])
	densubsample <- density(coeffs[,2])
	densdiff <- density(coeffs[,3])
	
	X11(width=10, height=8)
	par(mfrow=c(2,2))
	plot( densample , col = "green", main = "Coefficient values from the Cox model")
	lines( densubsample )
	legend("topleft",legend = c("sample","subsample"), col = c("black","green"), pch = 15)

	hist(coeffs[,3], main = "Sample minus Subsample Coefficients", xlab="Coef Diff", col = "lightblue")
	rug(coeffs[,3])
	par(new=TRUE)
	plot(densdiff$x, densdiff$y, type="l", axes = FALSE, ylab="",xlab="", col = "red" )
	axis(4)
	
	boxplot(coeffs[,3], horizontal = TRUE, col = "lightblue", main = " Difference of coefficients - Summary ")
	
	X11(width=10, height=8)
	par(mfrow=c(2,1))
	plot(mainsample$fit, col = c(1,3), main = "Mainsample")
	legend("topright",legend = c("sample","subsample"), col = c("black","green"), pch = 15)

	if(cens == 1){
		plot(mainsample$fit, col = c(1,3), main = "Mainsample and subsample chosen")
		lines(subsample$fit2, col = c(1,3),lty = 2)
		legend("topright",legend = c("T1","T2"), col = c("black","green"), pch = 15)
	}

	if(cens == 0){
		plot(mainsample$fitnocens, col = c(1,3))
		lines(subsample$fit2, col = c(1,3),lty = 2)
		legend("topright",legend = c("T1","T2"), col = c("black","green"), pch = 15)
	}

		### Summary

	summary(coeffs[,3])
	#l <- list(mainsample=mainsample, subsample=subsample, coeffs=coeffs)
	#return(l)
}

fonction subsampling pour l'example

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
subsampling <- function(p1, p2, cens){
 
	cens <- 0
 
	mainsample1 <- cbind(mainsample$y1, mainsample$status[1:mainsample$n1], mainsample$treat[1:mainsample$n1])
	colnames(mainsample1) <- c("time","status","treat")
	mainsample1 <- data.frame(mainsample1)
 
	mainsample2 <- cbind(mainsample$y2, mainsample$status[mainsample$n1+1:mainsample$n2], mainsample$treat[mainsample$n1+1:mainsample$n2])
	colnames(mainsample2) <- c("time","status","treat")
	mainsample2 <- data.frame(mainsample2)
 
	y1tild <- mainsample1[sample(1:mainsample$n1,ceiling(p1*mainsample$n1), replace=FALSE),]
	y2tild <- mainsample2[sample(1:mainsample$n2,ceiling(p2*mainsample$n2), replace=FALSE),]
 
	n1 <- length(y1tild )
	n2 <- length(y2tild )
 
	ytild <- rbind(y1tild, y2tild)
 
	fit2 <- survfit(Surv(time, status) ~ treat, data=ytild)
	#plot(fit2, col = c(1,3))
	#legend("topleft",legend = c("T1","T2"), col = c("black","green"), pch = 15)
 
	l <- list(y1tild=y1tild, y2tild=y2tild, fit2=fit2, time=ytild$time, status=ytild$status, treat=ytild$treat,n1=n1, n2=n2)
	return(l)
 
}

pour appeler la fonction :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
sim(nb = 500,N = 600,r = 0.7,p1 = 0.6,p2 = 0.6, subsamp = subsampling, cens = 0)

Faites moi savoir si je ne suis pas assez claire ou s'il y a besoin de précisions.
Merci d'avance

[andy18]

C'est bon problème résolu. Il fallait ajouter l'attribut "mainsample" à la fonctio"subsampling".