|
9 | 9 | #' @include mixed_effect_class.R HSD_class.R
|
10 | 10 | #' @export HSDEM
|
11 | 11 | HSDEM<-setClass(
|
12 |
| - "HSDEM", |
13 |
| - contains=c('method','stato'), |
14 |
| - slots=c( |
15 |
| - # INPUTS |
16 |
| - params.alpha='entity.stato', |
17 |
| - params.mtc='entity.stato', |
18 |
| - params.formula='entity', |
19 |
| - |
20 |
| - # OUTPUTS |
21 |
| - outputs.p_value='entity.stato', |
22 |
| - outputs.significant='entity' |
23 |
| - ), |
24 |
| - prototype = list(name='Tukey Honest Significant Difference using estimated marginal means', |
25 |
| - description='Tukey HSD post hoc tests for mixed effects models using estimated marginal means', |
26 |
| - type="univariate", |
27 |
| - predicted='p_value', |
28 |
| - stato.id="STATO:0000187", |
29 |
| - |
30 |
| - params.alpha=entity.stato(name='Confidence level', |
31 |
| - stato.id='STATO:0000053', |
32 |
| - value=0.05, |
33 |
| - type='numeric', |
34 |
| - description='the p-value cutoff for determining significance.' |
| 12 | + "HSDEM", |
| 13 | + contains=c('method','stato'), |
| 14 | + slots=c( |
| 15 | + # INPUTS |
| 16 | + params.alpha='entity.stato', |
| 17 | + params.mtc='entity.stato', |
| 18 | + params.formula='entity', |
| 19 | + |
| 20 | + # OUTPUTS |
| 21 | + outputs.p_value='entity.stato', |
| 22 | + outputs.significant='entity' |
35 | 23 | ),
|
36 |
| - params.mtc=entity.stato(name='Multiple Test Correction method', |
37 |
| - stato.id='OBI:0200089', |
38 |
| - value='none', |
39 |
| - type='numeric', |
40 |
| - description='The method used to adjust for multiple comparisons.' |
41 |
| - ), |
42 |
| - |
43 |
| - outputs.p_value=entity.stato(name='p value', |
44 |
| - stato.id='STATO:0000175', |
45 |
| - type='numeric', |
46 |
| - description='the probability of observing the calculated t-statistic.' |
47 |
| - ), |
48 |
| - outputs.significant=entity(name='Significant features', |
49 |
| - #stato.id='STATO:0000069', |
50 |
| - type='logical', |
51 |
| - description='TRUE if the calculated p-value is less than the supplied threhold (alpha)' |
| 24 | + prototype = list(name='Tukey Honest Significant Difference using estimated marginal means', |
| 25 | + description='Tukey HSD post hoc tests for mixed effects models using estimated marginal means', |
| 26 | + type="univariate", |
| 27 | + predicted='p_value', |
| 28 | + stato.id="STATO:0000187", |
| 29 | + |
| 30 | + params.alpha=entity.stato(name='Confidence level', |
| 31 | + stato.id='STATO:0000053', |
| 32 | + value=0.05, |
| 33 | + type='numeric', |
| 34 | + description='the p-value cutoff for determining significance.' |
| 35 | + ), |
| 36 | + params.mtc=entity.stato(name='Multiple Test Correction method', |
| 37 | + stato.id='OBI:0200089', |
| 38 | + value='none', |
| 39 | + type='numeric', |
| 40 | + description='The method used to adjust for multiple comparisons.' |
| 41 | + ), |
| 42 | + |
| 43 | + outputs.p_value=entity.stato(name='p value', |
| 44 | + stato.id='STATO:0000175', |
| 45 | + type='numeric', |
| 46 | + description='the probability of observing the calculated t-statistic.' |
| 47 | + ), |
| 48 | + outputs.significant=entity(name='Significant features', |
| 49 | + #stato.id='STATO:0000069', |
| 50 | + type='logical', |
| 51 | + description='TRUE if the calculated p-value is less than the supplied threhold (alpha)' |
| 52 | + ) |
52 | 53 | )
|
53 |
| - ) |
54 | 54 | )
|
55 | 55 |
|
56 | 56 | #' @export
|
57 | 57 | setMethod(f="method.apply",
|
58 |
| - signature=c("HSDEM",'dataset'), |
59 |
| - definition=function(M,D) { |
60 |
| - X=dataset.data(D) |
61 |
| - lmer_formula=aov2lme(M$formula) |
62 |
| - var_names=all.vars(M$formula) |
63 |
| - var_names_1=var_names[1] |
64 |
| - var_names=var_names[-1] |
65 |
| - y=dataset.sample_meta(D)[var_names] |
66 |
| - |
67 |
| - # set the contrasts |
68 |
| - O=options('contrasts') # keep the old ones |
69 |
| - options(contrasts = c("contr.sum","contr.poly")) |
70 |
| - |
71 |
| - # attempt to detect within factors |
72 |
| - within=which(var_names %in% all.names(M$formula)[which('Error'== all.names(M$formula))+2]) |
73 |
| - if (length(within)>0) { |
74 |
| - var_names_ex=var_names[-within] |
75 |
| - } else { |
76 |
| - var_names_ex=var_names |
77 |
| - } |
78 |
| - |
79 |
| - FF=full_fact(var_names_ex) |
80 |
| - FF=apply(FF,1,function(x) var_names_ex[x==1]) |
81 |
| - FF=FF[-1] |
82 |
| - |
83 |
| - output=apply(X,2,function(x) { |
84 |
| - temp=y |
85 |
| - temp[[var_names_1]]=scale(x,center = TRUE,scale=TRUE) |
86 |
| - |
87 |
| - dona=FALSE |
88 |
| - |
89 |
| - testlm=tryCatch({ # if any warnings/messages set p-values to NA as unreliable |
90 |
| - LM=lme(lmer_formula$f,random=lmer_formula$random,method='ML',data=temp,na.action=na.omit) |
91 |
| - }, warning = function(w) { |
92 |
| - NA |
93 |
| - }, message = function(m) { |
94 |
| - NA |
95 |
| - }, error = function(e) { |
96 |
| - NA |
97 |
| - }) |
98 |
| - |
99 |
| - output2=list() |
100 |
| - for (k in 1:length(FF)) { |
101 |
| - if (!is.na(testlm[[1]])) { |
102 |
| - testhsd=tryCatch({ |
103 |
| - output2[[k]]=as.data.frame(pairs(emmeans(LM,FF[[k]],data=temp))) |
104 |
| - }, warning = function(w) { |
105 |
| - NA |
106 |
| - } , message = function(m) { |
107 |
| - NA |
108 |
| - }, error = function(e) { |
109 |
| - NA |
110 |
| - }) |
111 |
| - |
112 |
| - if (!is.data.frame(testhsd[1])) { |
113 |
| - output2[[k]]=NA |
114 |
| - } |
| 58 | + signature=c("HSDEM",'dataset'), |
| 59 | + definition=function(M,D) { |
| 60 | + X=dataset.data(D) |
| 61 | + lmer_formula=aov2lme(M$formula) |
| 62 | + var_names=all.vars(M$formula) |
| 63 | + var_names_1=var_names[1] |
| 64 | + var_names=var_names[-1] |
| 65 | + y=dataset.sample_meta(D)[var_names] |
| 66 | + |
| 67 | + # set the contrasts |
| 68 | + O=options('contrasts') # keep the old ones |
| 69 | + options(contrasts = c("contr.sum","contr.poly")) |
| 70 | + |
| 71 | + # attempt to detect within factors |
| 72 | + within=which(var_names %in% all.names(M$formula)[which('Error'== all.names(M$formula))+2]) |
| 73 | + if (length(within)>0) { |
| 74 | + var_names_ex=var_names[-within] |
115 | 75 | } else {
|
116 |
| - output2[[k]]=NA |
| 76 | + var_names_ex=var_names |
117 | 77 | }
|
118 |
| - } |
119 |
| - |
120 |
| - return(output2) |
121 |
| - }) |
122 |
| - |
123 |
| - p_value=lapply(output,function(x) { |
124 |
| - x=as.data.frame(x) |
125 |
| - return(x$p.value) |
126 |
| - }) |
127 |
| - ln=length(p_value[[1]]) |
128 |
| - p_value=lapply(p_value,function(x){ |
129 |
| - if (length(x)!=ln) { |
130 |
| - return(p_value[[1]]*NA) |
131 |
| - } else { |
132 |
| - return(x) |
133 |
| - } |
134 |
| - }) |
135 |
| - |
136 |
| - p_value=do.call("rbind",p_value) |
137 |
| - p_value=data.frame(p_value) |
138 |
| - colnames(p_value)=as.data.frame(output[[1]])$contrast |
139 |
| - |
140 |
| - # fdr correct |
141 |
| - M$p_value=apply(p_value,2,p.adjust,method=M$mtc) |
142 |
| - |
143 |
| - M$significant=M$p_value<M$alpha |
144 |
| - |
145 |
| - # reset contrasts |
146 |
| - options(O) |
147 |
| - |
148 |
| - return(M) |
149 |
| - } |
| 78 | + |
| 79 | + FF=full_fact(var_names_ex) |
| 80 | + FF=apply(FF,1,function(x) var_names_ex[x==1]) |
| 81 | + FF=FF[-1] |
| 82 | + |
| 83 | + output=apply(X,2,function(x) { |
| 84 | + temp=y |
| 85 | + temp[[var_names_1]]=scale(x,center = TRUE,scale=TRUE) |
| 86 | + |
| 87 | + dona=FALSE |
| 88 | + |
| 89 | + testlm=tryCatch({ # if any warnings/messages set p-values to NA as unreliable |
| 90 | + LM=lme(lmer_formula$f,random=lmer_formula$random,method='ML',data=temp,na.action=na.omit) |
| 91 | + }, warning = function(w) { |
| 92 | + NA |
| 93 | + }, message = function(m) { |
| 94 | + NA |
| 95 | + }, error = function(e) { |
| 96 | + NA |
| 97 | + }) |
| 98 | + |
| 99 | + output2=list() |
| 100 | + for (k in 1:length(FF)) { |
| 101 | + if (!is.na(testlm[[1]])) { |
| 102 | + testhsd=tryCatch({ |
| 103 | + output2[[k]]=as.data.frame(pairs(emmeans(LM,FF[[k]],data=temp))) |
| 104 | + }, warning = function(w) { |
| 105 | + NA |
| 106 | + } , message = function(m) { |
| 107 | + NA |
| 108 | + }, error = function(e) { |
| 109 | + NA |
| 110 | + }) |
| 111 | + |
| 112 | + if (!is.data.frame(testhsd[1])) { |
| 113 | + output2[[k]]=NA |
| 114 | + } |
| 115 | + } else { |
| 116 | + output2[[k]]=NA |
| 117 | + } |
| 118 | + } |
| 119 | + |
| 120 | + return(output2) |
| 121 | + }) |
| 122 | + |
| 123 | + p_value=lapply(output,function(x) { |
| 124 | + x=as.data.frame(x) |
| 125 | + return(x$p.value) |
| 126 | + }) |
| 127 | + ln=length(p_value[[1]]) |
| 128 | + p_value=lapply(p_value,function(x){ |
| 129 | + if (length(x)!=ln) { |
| 130 | + return(p_value[[1]]*NA) |
| 131 | + } else { |
| 132 | + return(x) |
| 133 | + } |
| 134 | + }) |
| 135 | + |
| 136 | + p_value=do.call("rbind",p_value) |
| 137 | + p_value=data.frame(p_value) |
| 138 | + colnames(p_value)=as.data.frame(output[[1]])$contrast |
| 139 | + |
| 140 | + # fdr correct |
| 141 | + M$p_value=apply(p_value,2,p.adjust,method=M$mtc) |
| 142 | + |
| 143 | + M$significant=M$p_value<M$alpha |
| 144 | + |
| 145 | + # reset contrasts |
| 146 | + options(O) |
| 147 | + |
| 148 | + return(M) |
| 149 | + } |
150 | 150 | )
|
151 | 151 |
|
152 | 152 |
|
|
0 commit comments