-
Notifications
You must be signed in to change notification settings - Fork 5
/
testMitec.R
161 lines (142 loc) · 5.94 KB
/
testMitec.R
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
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#' TestM.ITEC
#'
#' This function performs TestM.ITEC
#' @param X is a matrix of encounter histories with K occasions
#' @param freq is a vector of the number of individuals with the corresponding encounter history
#' @param verbose controls the level of the details in the outputs; default is TRUE for all details
#' @param rounding is the level of rounding for outputs; default is 3
#' @return This function returns a list with first component the overall test and second component a data.frame with occasion, the value of the test statistic, degree of freedom, p-value and test performed (chi-square, Fisher or none).
#' @author Olivier Gimenez <olivier.gimenez@cefe.cnrs.fr>, Rémi Choquet, Roger Pradel
#' @keywords package
#' @export
#' @examples
#' \donttest{
#' # Read in Geese dataset:
#' geese = system.file("extdata", "geese.inp", package = "R2ucare")
#' geese = read_inp(geese)
#'
#' # Get encounter histories and number of individuals with corresponding histories
#' geese.hist = geese$encounter_histories
#' geese.freq = geese$sample_size
#'
#' # perform TestM.ITEC
#' testMitec(geese.hist,geese.freq)
#' }
testMitec <- function(X,freq,verbose=TRUE,rounding=3){
# test of multinomial mixture / MMLM method Yantis et al.
# adaptation of program Pascal mix
#derocc=1-min(filtre);
#filtre=find(filtre);
# various quantities to define
k = ncol(X)
his = X
a = max(his)
# initialization
table_multi_mitec = data.frame(occasion = rep(NA,k-3),stat = rep(NA,k-3), df = rep(NA,k-3), p_val = rep(NA,k-3), test_perf = rep(FALSE,k-3))
marray = multimarray(X,freq)
debutligne = 1
finligne = seq(a,a*k,by=a)
debutcolonne = seq(1,a*(k-1),by=a)
fincolonne = a * (k-1)
datat = marray[,2:(ncol(marray)-1)] # extrait du m-array avec les revus, sans les relaches ni les jamais revus
#for (i in 2:(k-2)+derocc){ # boucle sur les occasions
for (i in 2:(k-2)){ # boucle sur les occasions
mixandbases = datat[debutligne:finligne[i],debutcolonne[i]:fincolonne]
for (j in 1:(i-2)){
if ((i-2)<1) break
mixandbases[1:a,] = mixandbases[1:a,] + mixandbases[(a+1):(2*a),]
mixandbases = mixandbases[-((a+1):(2*a)),]
}
for (j in 1:((ncol(mixandbases)-a)/a-1)){
if (((ncol(mixandbases)-a)/a-1)<1) break
mixandbases[,(a+1):(2*a)] = mixandbases[,(a+1):(2*a)] + mixandbases[,(2*a+1):(3*a)]
mixandbases = mixandbases[,-((2*a+1):(3*a))]
}
nk = nrow(mixandbases)
tabtemp = mixandbases[(nk-a+1):nk,]
#tabtemp = tabtemp[filtre,] # les bases filtrees
mixandbases = rbind(mixandbases[1:(nk-a),],tabtemp)
nj = nrow(tabtemp)
nk = nrow(mixandbases)
if (any(apply(tabtemp,1,sum)==0) | (nj==0)){
table_multi_mitec[i-1,1] = i
table_multi_mitec[i-1,2] = 0
table_multi_mitec[i-1,3] = 0
table_multi_mitec[i-1,4] = 0
table_multi_mitec[i-1,5] = 'None'
break
}
# pooling
mixandbases = pooling_mixtures(nk,nj,a,mixandbases)
# DEBUT DU CALCUL DU TEST DE MELANGE (traduit de Yantis et al)
indic_mixbasis = c(matrix(0,nrow=nrow(mixandbases)-nj,ncol=1),rep(1,nj)) # indique 1 ou 0 suivant que base ou melange
data = cbind(mixandbases,indic_mixbasis)
#if verbosity>=3
# strtable=[ strtable {strcat('occasion :',num2str(i))} ];
# taille=size(data,2);
# strtable=[ strtable {'--------------- Seen again - Seen again later '} ];
# for kk=1:size(data,1)
# if data(kk,taille)==0
# strtable=[ strtable {strcat('When last released | ',num2str(data(kk,1:taille-1))) }];
# else
# strtable=[ strtable {strcat('Currently released | ',num2str(data(kk,1:taille-1))) }];
# end
# end
#end
nk = nrow(data)
r = ncol(data)
ni = r - 1
data = t(data)
nature = data[r,] # il s'agit de indic_mixanbasis'
data = data[-r,] # on la supprime !!!
tri = which(nature!=0) # coordonnees des bases
nj = length(tri) # nombre de bases
tri = c(tri,which(nature==0)) # ajout des coordonnees des melanges
M = data[,tri] # on renumerote bases et melanges
totk = apply(M,2,sum) # effectif des colonnes
CoorMelVide = which(totk[(nj+1):nk]==0)
if (!(length(CoorMelVide)==0)) M = M[,-CoorMelVide]
nk = ncol(M)
totk = apply(M,2,sum) # actualisation des effectifs des colonnes
# Si aucune base n'est vide & si melanges
if (nj!=nk){
# NEW definition des bases
Np = t(M[,1:nj])
# definition des melanges
Mp = t(M[,(nj+1):nk])
# calcul des coefficients du melanges
res = coef_mixtures(Mp,Np)
Q = res$P
P = res$PI
A = res$GAM
Q = rbind(P,Q)
# calcul des valeurs attendues
theoriques = matrix(rep(totk,ni),byrow=T,nrow=ni) * t(Q)
# calcul du nombre de degres de liberte
df = (nk-nj)*(ni-nj)
# test LR
tempchi2 = gof_test(1,c(M),c(theoriques))
table_multi_mitec[i-1,1] = i
table_multi_mitec[i-1,2] = tempchi2
table_multi_mitec[i-1,3] = df
table_multi_mitec[i-1,4] = 1 - stats::pchisq(tempchi2,df)
table_multi_mitec[i-1,5] = 'Chi-square'
} else {
table_multi_mitec[i-1,1] = i
table_multi_mitec[i-1,2] = 0
table_multi_mitec[i-1,3] = 0
table_multi_mitec[i-1,4] = 0
table_multi_mitec[i-1,5] = 'None'
}
}
# compute overall test:
stat = sum(as.numeric(table_multi_mitec[,2]))
stat = round(stat,rounding)
dof = sum(as.numeric(table_multi_mitec[,3]))
pval = 1 - stats::pchisq(stat,dof)
pval = round(pval,rounding)
# if user specifies all outputs
if (verbose==TRUE) return(list(testMitec=c(stat=stat,df=dof,p_val=pval),details=table_multi_mitec))
# otherwise
if (verbose==FALSE) return(list(testMitec=c(stat=stat,df=dof,p_val=pval)))
}