Mercurial > repos > stevecassidy > parseeval

comparison parseval.R @ 0:6c1e450b02c3 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 72cee9103c0ae4acb5794afaed179bea2c729f2c-dirty
author stevecassidy
date Sat, 11 Mar 2017 21:35:38 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6c1e450b02c3
1 parseEval <- function (c1, c2, c3, cipi12, cipi23, voweld, anchor=15, var_anchor=3, types, words=12, simN=1000, RE_rsd, CC_rsd, LE_rsd) {
2 # basic validy-check ... input parameters
3 c1<-as.numeric(c1)
4 c2<-as.numeric(c2)
5 c3<-as.numeric(c3)
6 cipi12<-as.numeric(cipi12)
7 cipi23<-as.numeric(cipi23)
8 voweld<-as.numeric(voweld)
9 anchor<-as.integer(anchor)
10 var_anchor<-as.integer(var_anchor)
11 types<-as.integer(types)
12 words<-as.integer(words)
13 simN<-as.integer(simN)
14 RE_rsd<-as.numeric(RE_rsd)
15 CC_rsd<-as.numeric(CC_rsd)
16 LE_rsd<-as.numeric(LE_rsd)
17 # phonetic parameters
18 C1p<-c3[1] # plateau duration of the first consonant in triconsonantal clusters
19 C1stdv<-c3[2] # standard deviation of the first consonant in triconsonantal clusters
20 C2p<-c2[1] # plateau duration of the second consonant in triconsonantal clusters and the first consonant in bi-consonantal clusters
21 C2stdv<-c2[2] # standard deviation of the second consonant in triconsonantal clusters and the first consonant in bi-consonantal clusters
22 C3p<-c1[1] # plateau duration of the immediately prevocalic consonant
23 C3stdv<-c1[2] # standard deviation of the immediately prevocalic consonant
24 C12ipi<-cipi12[1] # the duration of the interval between the plateaus of the first two consonants in triconsonantal clusters
25 C12stdv<-cipi12[2] # the standard deviation of the interval between the plateaus of the first two consonants in triconsonantal clusters
26 C23ipi<-cipi23[1] # the duration of the interval between the plateaus of the first two consonants in biconsonantal clusters and between the second two consonants in triconsonantal clusters
27 C23stdv<-cipi23[2] # the standard deviation of the interval between the plateaus of the first two consonants in biconsonantal clusters and between the second two consonants in triconsonantal clusters
28 vowel_duration<-voweld # the duration of the vowel
29 variability_range<-anchor # number of stepwise increases in variability simulated by the model
30 variability_resolution<-var_anchor # size of each stepwise increase in variability
31 words_per_word_type<-words # the number of words (stimuli) per word type, aka "little n"
32 word_types<-types # number of different word types, e.g. #CV-, #CCV-, #CCCV-
33 data_RSD<-c(RE_rsd, LE_rsd, CC_rsd) # lumps RSD measures into single array
34 #creating matrices for later use
35 A_simp <- matrix(nrow=variability_range, ncol=words_per_word_type)
36 A_comp <- matrix(nrow=variability_range, ncol=words_per_word_type)
37 # creating matrices to hold the SD values
38 LE_SD_simp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
39 LE_SD_comp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
40 RE_SD_simp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
41 RE_SD_comp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
42 CC_SD_simp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
43 CC_SD_comp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
44 # creating matrices to hold the RSD values
45 LE_RSD_simp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
46 LE_RSD_comp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
47 RE_RSD_simp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
48 RE_RSD_comp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
49 CC_RSD_simp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
50 CC_RSD_comp<-matrix(nrow=simN, ncol=variability_range, byrow=TRUE)
51 if (word_types==3) {
52 tepa<-c("Testing Triads")
53 # print(c("Simulating Data (Triads)"), quote=F)
54 # pb <- txtProgressBar(min = 0, max = simN, style = 3)
55 for (count in 1:simN) {
56 # setTxtProgressBar(pb, count)
57 # generate CCC tokens
58 # generate timestamps for C3 (the prevocalic consonant)
59 # generate general error-term for C3
60 e <- C3stdv*(rnorm(1))
61 # generate R(ight plateau edge = Release) of prevocalic consonant
62 # generate words_per_word_type/word_types Gaussian distributed numbers (for CCC tokens only)
63 # with mean 500, variance 10
64 CCCR3 <- rnorm(words_per_word_type/word_types, mean=500, sd=sqrt(20))
65 # generate L(eft plateau edge = Target) of prevocalic consonant
66 CCCL3 <- CCCR3 - C3p + e #generate L3 corresponding to R3 by assuming a plateau duration of C3p
67 # calculate midpoint of prevocalic consonant plateau
68 CCCM3 <- ((CCCR3 + CCCL3) / 2)
69 # generate timestamps for C2
70 # generate general error-term for C2
71 e1 <- C23stdv * (rnorm(1)) #normally distributed random error mean=0, sd=1
72 e2 <- C2stdv * (rnorm(1)) #normally distributed random error mean=0, sd=1
73 # Generate right edge of C2
74 CCCR2 <- CCCL3 - C23ipi + e1 # generate right edge of C2 from left edge of C3 assuming an ipi of C23ipi
75 # Generate left edge of C2
76 CCCL2 <- CCCR2 - C2p + e2 # generate left edge from right edge by assuming a plateau duration
77 # Calculate midpoint of C2
78 CCCM2 <- ((CCCR2+CCCL2)/2)
79 # generate timestamps for C1
80 # generate general error-term for C1
81 e1 <- C12stdv * (rnorm(1)) # normally distributed random error
82 e2 <- C3stdv * (rnorm(1))
83 # Generate right edge
84 CCCR1 <- CCCL2 - C12ipi + e1 # generate right edge of C1 from left edge of C2 assuming ipi of 40ms
85 # generate L(eft plateau edge = Target) of C1
86 CCCL1 <- CCCR1 - C3p + e2 # generate L2 corresponding to CR1 by assuming a plateau of 10ms
87 # calculate midpoint of prevocalic consonant
88 CCCM1 <- ((CCCR1 + CCCL1)/2) # right edge of C1
89 #generate CC tokens
90 #generate timestamps for C3 (prevocalic consonant)
91 # generate general error-term for C3
92 e <- C3stdv * (rnorm(1)) # normally distributed random error
93 # generate R(ight plateau edge = Release) of prevocalic consonant
94 CCR3 <- rnorm(words_per_word_type/word_types, mean=500, sd=sqrt(20)) # generate N Gaussian distributed numbers with mean 500, variance 10
95 # generate L(eft plateau edge = Target) of prevocalic consonant
96 CCL3 <- CCR3 - C3p + e # generate L3 corresponding to R3 by assuming a plateau duration of C3p
97 # calculate midpoint of prevocalic consonant plateau
98 CCM3 <- ((CCR3 + CCL3) / 2)
99 #generate timestamps for C2
100 # generate general error-term for C2
101 e1 <- C23stdv * (rnorm(1))
102 e2 <- C2stdv * (rnorm(1))
103 # Generate right edge of C2
104 CCR2 <- CCL3 - C23ipi + e1 # generate right edge of C2 from left edge of C3 assuming an ipi of C23ipi
105 # Generate left edge of C2
106 CCL2 <- CCR2 - C2p + e2 # generate left edge from right edge by assuming a plateau duration
107 # Calculate midpoint of C2
108 CCM2 <- ((CCR2 + CCL2) / 2)
109 # generate C tokens
110 # generate timestamps for C3 (the prevocalic consonant)
111 # generate general error-term for C3
112 e <- C3stdv * (rnorm(1))
113 # Generate R(ight plateau edge = Release) of prevocalic consonant
114 CR3 <- rnorm(words_per_word_type/word_types, mean=500, sd=sqrt(20)) # generate N Gaussian distributed numbers with mean 500, variance 10
115 # generate L(eft plateau edge = Target) of prevocalic consonant
116 CL3 <- CR3 - C3p + e # generate L3 corresponding to R3 by assuming a plateau duration of C3p
117 # calculate midpoint of prevocalic consonant plateau
118 CM3 <- ((CR3 + CL3) / 2)
119 # generate timestamps for CCglobal
120 # for CCC clusters
121 CCglobal <- apply(cbind(CCCM1, CCCM2, CCCM3), 1, mean) #mean of consonant plateaux midpoints
122 # for CC clusters
123 CCglobal <- append(CCglobal, apply(cbind(CCM2, CCM3), 1, mean)) # mean of consonant plateaux midpoints
124 # for C clusters
125 CCglobal <- append(CCglobal, CM3)
126 # populate a single array with the midpoint of the pre-vocalic
127 # consonant of every word type; this array will be used to generate anchors
128 # for CCC clusters
129 Global_CM3 <- CCCM3 # mean of consonant plateaux midpoints
130 # for CC clusters
131 Global_CM3 <- append(Global_CM3, CCM3) # mean of consonant plateaux midpoints
132 # for C clusters
133 Global_CM3 <- append(Global_CM3, CM3)
134 # populate a single array with the Left_edge of the consonant cluster for every token
135 # this array will be used to calculate SD and RSD for EDGE to Anchor intervals
136 # for CCC clusters
137 Global_CL1 <- CCCL1 # Assigns the left edge of tri-consonantal tokens to the first third of Global_Cl1
138 # for CC clusters
139 Global_CL1 <- append(Global_CL1, CCL2) # Assigns the left edge of bi-consonantal tokens to the second third of Global_Cl1
140 # for C clusters
141 Global_CL1 <- append(Global_CL1, CL3) # Assigns the left edge of mono-consonantal tokens to the last third of Global_Cl1
142 # populate a single array with the Right_edge of the consonant cluster for every token
143 # this array is used to calculate SD and RSD for EDGE to Anchor intervals
144 # for CCC clusters
145 Global_CR3 <- CCCR3 # mean of consonant plateaux midpoints
146 # for CC clusters
147 Global_CR3 <- append(Global_CR3, CCR3) # mean of consonant plateaux midpoints
148 # for C clusters
149 Global_CR3 <- append(Global_CR3, CR3) # CCglobal synchronous with prevocalic consonant's plateau midpoint
150 # generate series of anchor points increasing in variability and/or distance from
151 # the prevocalic consonant reset the anchor array to zero
152 # one row for each anchor and one column for each token
153
154 # loop produces data/anchor for each token based on Simplex Hypothesis
155 stdv <- 0 # reset the value of the anchor stdev to zero
156 Ae <- NULL # reset anchor-error-term
157 for (cycle in 1:variability_range){ #creates multiple anchor points for each token
158 for (m in 1:words_per_word_type){ #creates anchor point for each token from the right edge of the token
159 Ae<-stdv*(rnorm(n=1)) #normally distributed random error, assuming mean of 0
160 A_simp[cycle, m] <- Global_CM3[m] + vowel_duration + Ae #generate anchor A according to the simplex onset hypothesis
161 }
162 stdv<-stdv+variability_resolution #creates new anchor point
163 }
164 # loop produces data/anchor for each token based on Complex Hypothesis
165 stdv <- 0 # reset the value of the anchor stdev to zero
166 Ae <- NULL # reset anchor-error-term
167 for (cycle in 1:variability_range) { #creates multiple anchor points for each token
168 for (m in 1:words_per_word_type) { #creates anchor point for each token from the right edge of the token
169 Ae<-stdv*(rnorm(1)) #normally distributed random error, assuming mean of 0
170 A_comp[cycle, m]<-CCglobal[m]+vowel_duration+Ae #generate anchor A according to the complex onset hypothesis
171 }
172 stdv<-stdv+variability_resolution #creates new anchor point
173 }
174 # Note about consonantal landmarks:
175 # they are replaced with each cycle of the simulation
176 # in constrast, RSD values for each landmark are stored across simulations.
177 # creating matrices to hold the SD values
178 x <- function(x) {sd(x-Global_CL1)}
179 y <- function(y) {sd(y-Global_CR3)}
180 z <- function(z) {sd(z-CCglobal)}
181 # computing the SD values
182 LE_SD_simp[count,] <- apply(A_simp, 1, x)
183 LE_SD_comp[count,] <- apply(A_comp, 1, x)
184 RE_SD_simp[count,] <- apply(A_simp, 1, y)
185 RE_SD_comp[count,] <- apply(A_comp, 1, y)
186 CC_SD_simp[count,] <- apply(A_simp, 1, z)
187 CC_SD_comp[count,] <- apply(A_comp, 1, z)
188 # computing the RSD values
189 LE_RSD_simp[count,] <- (apply(A_simp, 1, x))/((apply(A_simp, 1, mean))-mean(Global_CL1))
190 LE_RSD_comp[count,] <- (apply(A_comp, 1, x))/((apply(A_comp, 1, mean))-mean(Global_CL1))
191 RE_RSD_simp[count,] <- (apply(A_simp, 1, y))/((apply(A_simp, 1, mean))-mean(Global_CR3))
192 RE_RSD_comp[count,] <- (apply(A_comp, 1, y))/((apply(A_comp, 1, mean))-mean(Global_CR3))
193 CC_RSD_simp[count,] <- (apply(A_simp, 1, z))/(apply(A_simp, 1, mean)-mean(CCglobal))
194 CC_RSD_comp[count,] <- (apply(A_comp, 1, z))/(apply(A_comp, 1, mean)-mean(CCglobal))
195 }
196 # close(pb)
197 }
198 if (word_types==2) {
199 tepa<-c("Testing Dyads")
200 # print(c("Simulating Data (Dyads)"), quote=F)
201 # pb <- txtProgressBar(min = 0, max = simN, style = 3)
202 for (count in 1:simN) {
203 # setTxtProgressBar(pb, count)
204 # generate CCC tokens
205 # generate timestamps for C3 (the prevocalic consonant)
206 # generate general error-term for C3
207 e <- C3stdv * (rnorm(1))
208 # generate R(ight plateau edge = Release) of prevocalic consonant
209 # generate words_per_word_type/word_types Gaussian distributed numbers (for CCC tokens only)
210 # with mean 500, variance 10
211 CCCR3 <- rnorm(words_per_word_type/word_types, mean=500, sd=sqrt(20))
212 # generate L(eft plateau edge = Target) of prevocalic consonant
213 CCCL3 <- abs(CCCR3 - C3p + e) #generate L3 corresponding to R3 by assuming a plateau duration of C3p
214 # calculate midpoint of prevocalic consonant plateau
215 CCCM3 <- abs((CCCR3 + CCCL3) / 2)
216 # generate timestamps for C2
217 # generate general error-term for C2
218 e1 <- C23stdv * (rnorm(1)) #normally distributed random error
219 e2 <- C2stdv * (rnorm(1)) #normally distributed random error
220 # Generate right edge of C2
221 CCCR2 <- abs(CCCL3 - C23ipi + e1) # generate right edge of C2 from left edge of C3 assuming an ipi of C23ipi
222 # Generate left edge of C2
223 CCCL2 <- abs(CCCR2 - C2p + e2) # generate left edge from right edge by assuming a plateau duration
224 # Calculate midpoint of C2
225 CCCM2 <- abs((CCCR2 + CCCL2) / 2)
226 # generate timestamps for C1
227 # generate general error-term for C1
228 e1 <- C12stdv * (rnorm(1)) # normally distributed random error
229 e2<-C3stdv * (rnorm(1))
230 # Generate right edge
231 CCCR1 <- abs(CCCL2 - C12ipi + e1) # generate right edge of C1 from left edge of C2 assuming ipi of 40ms
232 # generate L(eft plateau edge = Target) of C1
233 CCCL1 <- abs(CCCR1 - C3p + e2) # generate L2 corresponding to CR1 by assuming a plateau of 10ms
234 # calculate midpoint of prevocalic consonant
235 CCCM1 <- abs((CCCR1 + CCCL1) / 2) # right edge of C1
236 #generate CC tokens
237 #generate timestamps for C3 (prevocalic consonant)
238 # generate general error-term for C3
239 e <- C3stdv * (rnorm(1)) # normally distributed random error, 0 mean
240 # generate R(ight plateau edge = Release) of prevocalic consonant
241 CCR3 <- rnorm(words_per_word_type/word_types, mean=500, sd=sqrt(20)) # generate N Gaussian distributed numbers with mean 500, variance 10
242 # generate L(eft plateau edge = Target) of prevocalic consonant
243 CCL3 <- abs(CCR3 - C3p + e) # generate L3 corresponding to R3 by assuming a plateau duration of C3p
244 # calculate midpoint of prevocalic consonant plateau
245 CCM3 <- abs((CCR3 + CCL3) / 2)
246 #generate timestamps for C2
247 # generate general error-term for C2
248 e1 <- C23stdv * (rnorm(1))
249 e2 <- C2stdv * (rnorm(1))
250 # Generate right edge of C2
251 CCR2 <- abs(CCL3 - C23ipi + e) # generate right edge of C2 from left edge of C3 assuming an ipi of C23ipi
252 # Generate left edge of C2
253 CCL2 <- abs(CCR2 - C2p + e) # generate left edge from right edge by assuming a plateau duration
254 # Calculate midpoint of C2
255 CCM2 <- abs((CCR2 + CCL2) / 2)
256 # generate timestamps for CCglobal
257 # for CCC clusters
258 CCglobal <- apply(cbind(CCCM1, CCCM2, CCCM3), 1, mean) #mean of consonant plateaux midpoints
259 # for CC clusters
260 CCglobal <- append(CCglobal, apply(cbind(CCM2, CCM3), 1, mean)) # mean of consonant plateaux midpoints
261 # populate a single array with the midpoint of the pre-vocalic
262 # consonant of every word type; this array will be used to generate anchors
263 # for CCC clusters
264 Global_CM3 <- CCCM3 # mean of consonant plateaux midpoints
265 # for CC clusters
266 Global_CM3 <- append(Global_CM3, CCM3, after=length(CCCM3)) # mean of consonant plateaux midpoints
267 # populate a single array with the Left_edge of the consonant cluster for every token
268 # this array will be used to calculate SD and RSD for EDGE to Anchor intervals
269 # for CCC clusters
270 Global_CL1 <- CCCL1 # Assigns the left edge of tri-consonantal tokens to the first third of Global_Cl1
271 # for CC clusters
272 Global_CL1 <- append(Global_CL1, CCL2, after=length(CCCL1)) # Assigns the left edge of bi-consonantal tokens to the second third of Global_Cl1
273 # populate a single array with the Right_edge of the consonant cluster for every token
274 # this array is used to calculate SD and RSD for EDGE to Anchor intervals
275 # for CCC clusters
276 Global_CR3 <- CCCR3 # mean of consonant plateaux midpoints
277 # for CC clusters
278 Global_CR3 <- append(Global_CR3, CCR3, after=length(CCCR3)) # mean of consonant plateaux midpoints
279 # generate series of anchor points increasing
280 # in variability and/or distance from the prevocalic consonant
281 stdv <- 0
282 Ae <- NULL
283 for (cycle in 1:variability_range){ #creates multiple anchor points for each token
284 for (m in 1:words_per_word_type){ #creates anchor point for each token from the right edge of the token
285 Ae<-stdv*(rnorm(n=1)) #normally distributed random error, assuming mean of 0
286 A_simp[cycle, m]<-Global_CM3[m] + vowel_duration + Ae #generate anchor A according to the simplex onset hypothesis
287 }
288 stdv<-stdv+variability_resolution
289 }
290 # loop produces anchor for each token based on Complex Hypothesis
291 stdv <- 0
292 Ae <- NULL
293 for (cycle in 1:variability_range){ #creates multiple anchor points for each token
294 for (m in 1:words_per_word_type){ #creates anchor point for each token from the right edge of the token
295 Ae<-stdv*(rnorm(1)) #normally distributed random error, assuming mean of 0
296 A_comp[cycle, m]<-CCglobal[m]+vowel_duration+Ae #generate anchor A according to the complex onset hypothesis
297 }
298 stdv<-stdv+variability_resolution #creates new anchor point
299 }
300 #creating matrices to hold the SD values
301 x <- function(x) {sd(x-Global_CL1)}
302 y <- function(y) {sd(y-Global_CR3)}
303 z <- function(z) {sd(z-CCglobal)}
304 # computing the SD values
305 LE_SD_simp[count,] <- abs(apply(A_simp, 1, x))
306 LE_SD_comp[count,] <- abs(apply(A_comp, 1, x))
307 RE_SD_simp[count,] <- abs(apply(A_simp, 1, y))
308 RE_SD_comp[count,] <- abs(apply(A_comp, 1, y))
309 CC_SD_simp[count,] <- abs(apply(A_simp, 1, z))
310 CC_SD_comp[count,] <- abs(apply(A_comp, 1, z))
311 # computing the RSD values
312 LE_RSD_simp[count,] <- abs((apply(A_simp, 1, x))/((apply(A_simp, 1, mean))-mean(Global_CL1)))
313 LE_RSD_comp[count,] <- abs((apply(A_comp, 1, x))/((apply(A_comp, 1, mean))-mean(Global_CL1)))
314 RE_RSD_simp[count,] <- abs((apply(A_simp, 1, y))/((apply(A_simp, 1, mean))-mean(Global_CR3)))
315 RE_RSD_comp[count,] <- abs((apply(A_comp, 1, y))/((apply(A_comp, 1, mean))-mean(Global_CR3)))
316 CC_RSD_simp[count,] <- abs((apply(A_simp, 1, z))/(apply(A_simp, 1, mean)-mean(CCglobal)))
317 CC_RSD_comp[count,] <- abs((apply(A_comp, 1, z))/(apply(A_comp, 1, mean)-mean(CCglobal)))
318 }
319 }
320 # close(pb)
321 # pb <- txtProgressBar(min = 1, max = variability_range, style = 3)
322 # assorted variables for diagnostics / plotting
323 aip_1<-rep(c(1:variability_range), 3)
324 edgep_1<-rep(c("LE_RSD", "RE_RSD", "CC_RSD"), each=variability_range)
325 LE_RSD_simp_median<-apply(apply(LE_RSD_simp, 2, sort), 2, median)
326 RE_RSD_simp_median<-apply(apply(RE_RSD_simp, 2, sort), 2, median)
327 CC_RSD_simp_median<-apply(apply(CC_RSD_simp, 2, sort), 2, median)
328 LE_RSD_comp_median<-apply(apply(LE_RSD_comp, 2, sort), 2, median)
329 RE_RSD_comp_median<-apply(apply(RE_RSD_comp, 2, sort), 2, median)
330 CC_RSD_comp_median<-apply(apply(CC_RSD_comp, 2, sort), 2, median)
331 simp<-c(LE_RSD_simp_median, RE_RSD_simp_median, CC_RSD_simp_median)
332 comp<-c(LE_RSD_comp_median, RE_RSD_comp_median, CC_RSD_comp_median)
333 RE_RSD_median<-c(RE_RSD_simp_median, RE_RSD_comp_median)
334 CC_RSD_median<-c(CC_RSD_simp_median, CC_RSD_comp_median)
335 # median RDSs across simulations as a function of anchorindex
336 plot.1<-data.frame(anchorindex=aip_1, edge=edgep_1, parse_s=simp, parse_c=comp)
337 # aggregating data for goodness of fit evaluation
338 RE_RSD_simp<-t(RE_RSD_simp)
339 LE_RSD_simp<-t(LE_RSD_simp)
340 CC_RSD_simp<-t(CC_RSD_simp)
341 RE_RSD_comp<-t(RE_RSD_comp)
342 LE_RSD_comp<-t(LE_RSD_comp)
343 CC_RSD_comp<-t(CC_RSD_comp)
344 # looping through the data to get the gof results
345 data_simp<-matrix(ncol=4)
346 data_comp<-matrix(ncol=4)
347 sigfit<-function(x) {
348 if(x > 98.503) {
349 SigFit<-1
350 } else {
351 SigFit<-0
352 }
353 }
354 # analyzing data simplex
355 # print(c("Analysing... simple parse"), quote=F)
356 for (i in 1 : variability_range) {
357 # setTxtProgressBar(pb, i)
358 sim_RSD<-cbind(RE_RSD_simp[i,], LE_RSD_simp[i,],CC_RSD_simp[i,])
359 temp<-apply(sim_RSD, 1, function(x) (lm(data_RSD ~ x)))
360 # organizing data for final analyses
361 # creating anchor-index
362 anchor_idx<-rep(i, times=simN)
363 # extracting F-Statistics
364 fstat<-unlist(lapply(temp, function(x) summary(x)$fstatistic[1]))
365 # extracting R-Squared values
366 rsquared<-unlist(lapply(temp, function(x) summary(x)$r.squared))
367 # check for SigFit
368 sgf<-sapply(fstat, sigfit)
369 # aggregating data
370 agg_mat<-matrix(data=c(anchor_idx, fstat, rsquared, sgf), nrow=length(anchor_idx), ncol=4, dimnames=list(c(NULL), c("Anchorindex", "Fratio", "Rsquared", "SigFit")))
371 # adding sgf to existing data
372 data_simp<-rbind(data_simp, agg_mat)
373 }
374 outp_sp<-temp
375 data_simp<-data_simp[complete.cases(data_simp),]
376 data_simp<-as.data.frame(data_simp)
377 data_simp$Anchorindex<-as.factor(data_simp$Anchorindex)
378 output_simp<-tapply(data_simp$SigFit, data_simp$Anchorindex, sum)
379 # analyzing data complex
380 # print(c("Analysing... complex parse"), quote=F)
381 for (i in 1 : variability_range) {
382 # setTxtProgressBar(pb, i)
383 sim_RSD<-cbind(RE_RSD_comp[i,], LE_RSD_comp[i,],CC_RSD_comp[i,])
384 temp<-apply(sim_RSD, 1, function(x) (lm(data_RSD ~ x)))
385 # organizing data for final analyses
386 anchor_idx<-rep(i, times=simN)
387 # extracting F-Statistics
388 fstat<-unlist(lapply(temp, function(x) summary(x)$fstatistic[1]))
389 # extracting R-Squared values
390 rsquared<-unlist(lapply(temp, function(x) summary(x)$r.squared))
391 # check for SigFit
392 sgf<-sapply(fstat, sigfit)
393 # aggregating data
394 agg_mat<-matrix(data=c(anchor_idx, fstat, rsquared, sgf), nrow=length(anchor_idx), ncol=4, dimnames=list(c(NULL), c("Anchorindex", "Fratio", "Rsquared", "SigFit")))
395 # adding sgf to existing data
396 data_comp<-rbind(data_comp, agg_mat)
397 }
398 outp_cp<-temp
399 data_comp<-data_comp[complete.cases(data_comp),]
400 data_comp<-as.data.frame(data_comp)
401 data_comp$Anchorindex<-as.factor(data_comp$Anchorindex)
402 output_comp<-tapply(data_comp$SigFit, data_comp$Anchorindex, sum)
403 # diagnostic plot 2
404 output_plot.2<-cbind(output_simp, output_comp)
405 names(output_plot.2)<-NULL
406 colnames(output_plot.2)<-c("parse_s", "parse_c")
407 aip_2<-(1:variability_range)
408 plot.2<-data.frame(anchorindex=aip_2, output_plot.2, hitr_s=(output_simp/simN), hitr_c=(output_comp/simN))
409 # assessing overall model quality
410 # sum of hits per number of simulations
411 modq_s<-(sum(plot.2[,2]))/simN
412 modq_c<-(sum(plot.2[,3]))/simN
413 # assorted data for third diagnostic plot
414 # sorting by Rsquared (asc), tie-breaker by Fratio (asc)
415 data_simp_o<-data_simp[order(data_simp[,3], data_simp[,2]),]
416 data_comp_o<-data_comp[order(data_comp[,3], data_comp[,2]),]
417 aip_3<-rep(c(1:variability_range), 2)
418 parse.f<-rep(c("simp","comp"), each=variability_range)
419 # median
420 simp_rs_median<-tapply(data_simp_o$Rsquared, data_simp_o$Anchorindex, median)
421 comp_rs_median<-tapply(data_comp_o$Rsquared, data_comp_o$Anchorindex, median)
422 simp_fr_median<-tapply(data_simp_o$Fratio, data_simp_o$Anchorindex, median)
423 comp_fr_median<-tapply(data_comp_o$Fratio, data_comp_o$Anchorindex, median)
424 rs_median<-c(simp_rs_median, comp_rs_median)
425 fr_median<-c(simp_fr_median, comp_fr_median)
426 plot.3_median<-data.frame(anchorindex=aip_3, parse=parse.f, rs_median=rs_median, fr_median=fr_median)
427 # mean
428 simp_rs_mean<-tapply(data_simp_o$Rsquared, data_simp_o$Anchorindex, mean)
429 comp_rs_mean<-tapply(data_comp_o$Rsquared, data_comp_o$Anchorindex, mean)
430 simp_fr_mean<-tapply(data_simp_o$Fratio, data_simp_o$Anchorindex, mean)
431 comp_fr_mean<-tapply(data_comp_o$Fratio, data_comp_o$Anchorindex, mean)
432 rs_mean<-c(simp_rs_mean, comp_rs_mean)
433 fr_mean<-c(simp_fr_mean, comp_fr_mean)
434 plot.3_mean<-data.frame(anchorindex=aip_3, parse=parse.f, rs_mean=rs_mean, fr_mean=fr_mean)
435 # prepare for output
436 # close(pb)
437 output<-list("perf"=list("Performance Simple"=modq_s, "Performance Complex"=modq_c), "mode"=tepa, "Plot_1"=plot.1, "Plot_2"=plot.2, "Plot_3"=list("mean"=plot.3_mean, "median"=plot.3_median), "reg"=list("simplex_parse"<-outp_sp, "complex_parse"<-outp_cp), "sim_RSD"=list("simp"=data_simp, "comp"=data_comp))
438 cat("\n", "\n","Overall Quality of Modell-Performance", "\t", "(", tepa, ")", "\n",
439 "(Ratio of:","\t", "Total Number of Hits / Number of Simulations)","\n",
440 "------------------------","\n",
441 "Simple Modelling:", "\t", modq_s, "\t","\t","\t","\t", sum(plot.2[,2])," / ", simN, "\n", "\n",
442 "Complex Modelling:", "\t", modq_c, "\t","\t","\t","\t", sum(plot.2[,3])," / ", simN, "\n", "\n", sep="")
443 return(invisible(output))
444 }