Messieurs-Dames bonjour,
Je débute en R, essayant de comprendre du mieux que je peux l'utilisation de ses diverses fonctions, mais il y a certaines choses qui restent obscurs, c'est pourquoi je me tourne vers vous dans l'espoir que vous puissiez éclairer ma lanterne.
Je possède un jeu de données qui représente plusieurs évènements, avec, pour chacun d'eux, des coordonnées (X ; Y) ainsi qu'une date (JJ/MM/AAAA). Je cherche à obtenir un tableau de contingence m'indiquant le nombre de ces évènement qui ont eu lieu dans une certaine fourchette spatiale et temporelle de la forme:
1-2 jrs | 3-4 jrs | 5-6 jrs | ...
0M|
1-100M|
100-200M|
(Cela s’appelle un test de Knox, qui a été utilisé dans le but de mettre en évidence les groupement ou "clusters" de maladie, et sa dispersion, dans le cas d"une épidémie)
(ce test met, par ailleurs, en relation les valeurs obtenu avec des valeurs estimés (issu d'une itération de Monte Carlo par exemple) pour donner la significativité de chacune des valeurs...)
Ayant cherché longuement sur la toile je suis parvenu à trouver plusieurs fonctions de R permettant d'effectuer cette relation spatio-temporelle.
Néanmoins je ne parviens pas à les appliquer:
1 - http://www.niph.go.jp/soshiki/gijuts...ad/Rfunctions/
2 - http://finzi.psych.upenn.edu/R/libra...html/knox.html =>(fonction) https://r-forge.r-project.org/scm/vi...e&pathrev=1292
Dans cette première fonction lorsque je rentre mes données (x,y,time,del1,del2) et que j'éxecute la fonction, la console m'indique:
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86 KnoxA.test<-function(x,y,time,del1,del2){ # ----------------------------------------------------------------- # An example of execution # # dat<-scan("d:/demo/McHardy.dat",list(x=0,y=0,time=0)) # del1<-2; del2<-0 # out<-KnoxA.test(dat$x,dat$y,dat$time,del1,del2) # # --------------------------------------------------------------------- # ARGUEMENTS # x : a vector of x-coordinates of case # y : a vector of y-coordinates of case # time : a vector of observed times # del1 : a measure of closeness in space # del2 : a measure of closeness in time # # VALUES # Knox.T : test statistic # Expected : expected number of cases # Var : variance # Standardized : standardized test statistic # Poisson.p.value : p-value by Poisson approximation # Poisson.Mid.p.value : mid-p-value by Poisson approximation # Normal.p.value : p-value by Normal approximation # #----------------------------------------------------------------------- n<-length(x) # sdis<-matrix(0,n,n) tdis<-matrix(0,n,n) as<-sdis at<-tdis for (i in 1:n){ for (j in 1:n){ sdis[i,j]<-sqrt( (x[i]-x[j])^2 + (y[i]-y[j])^2 ) tdis[i,j]<- abs( time[i] - time[j]) } } # for(i in 1:n){ for (j in 1:n){ as[i,j] <- ifelse(sdis[i,j]<= del1, 1, 0) at[i,j] <- ifelse(tdis[i,j]<= del2, 1, 0) } as[i,i]<-0 at[i,i]<-0 } # s1<-0; s2<-0; s3<-0 for(i in 1:n){ for (j in 1:n){ s1<-s1+as[i,j] * at[i,j] s2<-s2+as[i,j] s3<-s3+at[i,j] } } obsT <- s1/2 Ms <- s2/2 Mt <- s3/2 # s1<-0; s2<-0 for(i in 1:n){ for (j in 1:n){ for (k in 1:n){ if(k!=j){ s1<-s1+as[i,j] * as[i,k] s2<-s2+at[i,j] * at[i,k] } } } } # Ms2 <- s1/2 Mt2 <- s2/2 ee <- Ms*Mt/(n*(n-1)/2) vv <- ee - ee*ee + 4*Ms*Mt/n/(n-1)/(n-2) vv <- vv + 4*(Ms*(Ms-1)-Ms2) * (Mt*(Mt-1)-Mt2)/n/(n-1)/(n-2)/(n-3) sdk <- sqrt(vv) zz <-(obsT - ee)/sdk pk <- 1-ppois(obsT-1,ee) mid <- 1-ppois(obsT,ee)+dpois(obsT,ee)/2 pkn <- 1-pnorm(zz) list(Knox.T=obsT, Expected=ee, Var=vv, Standardized=zz, Poisson.p.value=pk, Poisson.Mid.p.value=mid, Normal.p.value=pkn) }
Qu'est ce qui ne va pas?
Code : Sélectionner tout - Visualiser dans une fenêtre à part Error: unexpected '}' in "}"
Aussi, une deuxième fonction m'est proposé :
Mais pour cette fonction ci, je n'ai aucune idée de ce que je dois rentrer pour "dt,ds", ou plutôt, si j'ai bien compris, comment transformer mes X et Y pour le rentrer dans 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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138 ################################################################################ ### Part of the surveillance package, http://surveillance.r-forge.r-project.org ### Free software under the terms of the GNU General Public License, version 2, ### a copy of which is available at http://www.r-project.org/Licenses/. ### ### Knox test for space-time interaction ### ### Copyright (C) 2015 Sebastian Meyer ### $Revision$ ### $Date$ ################################################################################ knox <- function (dt, ds, eps.t, eps.s, simulate.p.value = TRUE, B = 999, ...) { stopifnot(length(dt) == length(ds)) ## logical vectors indicating which pairs are close in time and space closeInTime <- if (is.logical(dt)) { dt } else { stopifnot(is.numeric(dt), isScalar(eps.t)) dt <= eps.t } closeInSpace <- if (is.logical(ds)) { ds } else { stopifnot(is.numeric(ds), isScalar(eps.s)) ds <= eps.s } ## manually build the contingency table (table() with factor() is too slow) .lab <- c("close", "not close") knoxtab <- array( tabulate(4L - closeInTime - 2L*closeInSpace, nbins = 4L), dim = c(2L, 2L), dimnames = list( dt = if (is.logical(dt)) .lab else paste(c("<=", " >"), eps.t), ds = if (is.logical(ds)) .lab else paste(c("<=", " >"), eps.s) )) class(knoxtab) <- "table" ## expected number of close pairs in the absence of spatio-temporal interaction npairs <- sum(knoxtab) expected <- sum(knoxtab[1L,]) / npairs * sum(knoxtab[,1L]) ##<- this order of terms avoids integer overflow ## test statistic is the number of spatio-temporally close pairs METHOD <- "Knox test" STATISTIC <- knoxtab[1L] ## determine statistical significance pval_Poisson <- ppois(STATISTIC, expected, lower.tail = FALSE) PVAL <- if (simulate.p.value) { # Monte Carlo permutation approach stopifnot(isScalar(B)) B <- as.integer(B) METHOD <- paste(METHOD, "with simulated p-value") PARAMETER <- setNames(B, "B") permstats <- plapply(X = integer(B), FUN = function (...) sum(closeInSpace & closeInTime[sample.int(npairs)]), ...) structure(mean(c(STATISTIC, permstats, recursive = TRUE) >= STATISTIC), Poisson = pval_Poisson) } else { METHOD <- paste(METHOD, "with Poisson approximation") PARAMETER <- setNames(expected, "lambda") pval_Poisson } ## return test results structure( list(method = METHOD, data.name = paste("dt =", deparse(substitute(dt)), "and ds =", deparse(substitute(ds))), statistic = setNames(STATISTIC, "number of close pairs"), parameter = PARAMETER, p.value = PVAL, alternative = "greater", null.value = setNames(expected, "number"), permstats = if (simulate.p.value) { unlist(permstats, recursive = FALSE, use.names = FALSE) }, table = knoxtab), class = c("knox", "htest") ) } print.knox <- function (x, ...) { ## first print by the default method for class "htest" NextMethod("print") ## then also output the contingency table cat("contingency table:\n") print(x$table) cat("\n") invisible(x) } plot.knox <- function (x, ...) { if (is.null(permstats <- x[["permstats"]])) { stop("this plot-method is for a permutation-based Knox test") } defaultArgs <- list( permstats = permstats, xmarks = setNames(c(x[["null.value"]], x[["statistic"]]), c("expected", "observed")), xlab = "number of close pairs" ) do.call("epitestplot", modifyList(defaultArgs, list(...))) } xtable.knox <- function (x, caption = NULL, label = NULL, align = paste0("r|rr", if (!is.null(sumlabel)) "|r"), digits = 0, display = NULL, ..., sumlabel = "$\\sum$") { tab <- x$table if (!is.null(sumlabel)) { FUN <- setNames(list(sum), sumlabel) tab <- addmargins(tab, FUN = FUN, quiet = TRUE) } xtable(tab, caption = caption, label = label, align = align, digits = digits, display = display, ...) } toLatex.knox <- function (object, hline.after = NULL, sanitize.text.function = NULL, ...) { xtab <- xtable(object, ...) if (is.null(hline.after)) hline.after <- unique(c(-1,0,2,nrow(xtab))) if (is.null(sanitize.text.function)) sanitize.text.function <- function (x) gsub("<=", "$\\le$", gsub(">", "$>$", x, fixed = TRUE), fixed = TRUE) res <- print(xtab, hline.after = hline.after, sanitize.text.function = sanitize.text.function, ...) class(res) <- "Latex" invisible(res) # print.xtable() by default does 'print.results' }
(les détails de cette fonction sont indiqués ici => http://finzi.psych.upenn.edu/R/libra...html/knox.html).
ma question est donc la suivante:
Laquelle de ces deux fonctions vous semble la plus pertinente au vu de ce que je cherche à obtenir? et si aucune des deux le semble, quelles fonctions vous me conseilleriez d'explorer ?
Bien sur je ne cherche pas à ce que l'on effectue le travail à ma place, mais bien qu'ayant tant bien que mal cherché à me former en R, il y a certaines choses que je n'ai pas encore saisi et j'aurais besoin de votre avis ainsi que de vos conseils.
d'avance, merci pour l'intérêt que vous portez à ce message et pour votre réponse.
Cordialement,
JANIN Oscar
Partager