Mercurial > repos > dereeper > mlmm
comparison source_library/emma.r @ 1:380b364980f9 draft default tip
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
author | dereeper |
---|---|
date | Mon, 16 Apr 2018 08:50:05 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
0:6b7107812931 | 1:380b364980f9 |
---|---|
1 emma.kinship <- function(snps, method="additive", use="all") { | |
2 n0 <- sum(snps==0,na.rm=TRUE) | |
3 nh <- sum(snps==0.5,na.rm=TRUE) | |
4 n1 <- sum(snps==1,na.rm=TRUE) | |
5 nNA <- sum(is.na(snps)) | |
6 | |
7 stopifnot(n0+nh+n1+nNA == length(snps)) | |
8 | |
9 if ( method == "dominant" ) { | |
10 flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps)) | |
11 snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)] | |
12 } | |
13 else if ( method == "recessive" ) { | |
14 flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps)) | |
15 snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)] | |
16 } | |
17 else if ( ( method == "additive" ) && ( nh > 0 ) ) { | |
18 dsnps <- snps | |
19 rsnps <- snps | |
20 flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps)) | |
21 dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)] | |
22 flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps)) | |
23 rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)] | |
24 snps <- rbind(dsnps,rsnps) | |
25 } | |
26 | |
27 if ( use == "all" ) { | |
28 mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps)) | |
29 snps[is.na(snps)] <- mafs[is.na(snps)] | |
30 } | |
31 else if ( use == "complete.obs" ) { | |
32 snps <- snps[rowSums(is.na(snps))==0,] | |
33 } | |
34 | |
35 n <- ncol(snps) | |
36 K <- matrix(nrow=n,ncol=n) | |
37 diag(K) <- 1 | |
38 | |
39 for(i in 2:n) { | |
40 for(j in 1:(i-1)) { | |
41 x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j]) | |
42 K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x)) | |
43 K[j,i] <- K[i,j] | |
44 } | |
45 } | |
46 return(K) | |
47 } | |
48 | |
49 emma.eigen.L <- function(Z,K,complete=TRUE) { | |
50 if ( is.null(Z) ) { | |
51 return(emma.eigen.L.wo.Z(K)) | |
52 } | |
53 else { | |
54 return(emma.eigen.L.w.Z(Z,K,complete)) | |
55 } | |
56 } | |
57 | |
58 emma.eigen.L.wo.Z <- function(K) { | |
59 eig <- eigen(K,symmetric=TRUE) | |
60 return(list(values=eig$values,vectors=eig$vectors)) | |
61 } | |
62 | |
63 emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) { | |
64 if ( complete == FALSE ) { | |
65 vids <- colSums(Z)>0 | |
66 Z <- Z[,vids] | |
67 K <- K[vids,vids] | |
68 } | |
69 eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE) | |
70 return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE))) | |
71 } | |
72 | |
73 emma.eigen.R <- function(Z,K,X,complete=TRUE) { | |
74 if ( ncol(X) == 0 ) { | |
75 return(emma.eigen.L(Z,K)) | |
76 } | |
77 else if ( is.null(Z) ) { | |
78 return(emma.eigen.R.wo.Z(K,X)) | |
79 } | |
80 else { | |
81 return(emma.eigen.R.w.Z(Z,K,X,complete)) | |
82 } | |
83 } | |
84 | |
85 emma.eigen.R.wo.Z <- function(K, X) { | |
86 n <- nrow(X) | |
87 q <- ncol(X) | |
88 S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X) | |
89 eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE) | |
90 stopifnot(!is.complex(eig$values)) | |
91 return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)])) | |
92 } | |
93 | |
94 emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) { | |
95 if ( complete == FALSE ) { | |
96 vids <- colSums(Z) > 0 | |
97 Z <- Z[,vids] | |
98 K <- K[vids,vids] | |
99 } | |
100 n <- nrow(Z) | |
101 t <- ncol(Z) | |
102 q <- ncol(X) | |
103 | |
104 SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z) | |
105 eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE) | |
106 if ( is.complex(eig$values) ) { | |
107 eig$values <- Re(eig$values) | |
108 eig$vectors <- Re(eig$vectors) | |
109 } | |
110 qr.X <- qr.Q(qr(X)) | |
111 return(list(values=eig$values[1:(t-q)], | |
112 vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)), | |
113 complete=TRUE)[,c(1:(t-q),(t+1):n)])) | |
114 } | |
115 | |
116 emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) { | |
117 n <- length(xi) | |
118 delta <- exp(logdelta) | |
119 return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) ) | |
120 } | |
121 | |
122 emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) { | |
123 t <- length(xi.1) | |
124 delta <- exp(logdelta) | |
125 # stopifnot(length(lambda) == length(etas.1)) | |
126 return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) ) | |
127 } | |
128 | |
129 emma.delta.ML.dLL.wo.Z <- function(logdelta, lambda, etas, xi) { | |
130 n <- length(xi) | |
131 delta <- exp(logdelta) | |
132 etasq <- etas*etas | |
133 ldelta <- lambda+delta | |
134 return( 0.5*(n*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/(xi+delta))) ) | |
135 } | |
136 | |
137 emma.delta.ML.dLL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) { | |
138 t <- length(xi.1) | |
139 delta <- exp(logdelta) | |
140 etasq <- etas.1*etas.1 | |
141 ldelta <- lambda+delta | |
142 return( 0.5*(n*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/(xi.1+delta))+(n-t)/delta) ) ) | |
143 } | |
144 | |
145 emma.delta.REML.LL.wo.Z <- function(logdelta, lambda, etas) { | |
146 nq <- length(etas) | |
147 delta <- exp(logdelta) | |
148 return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas*etas/(lambda+delta))))-sum(log(lambda+delta))) ) | |
149 } | |
150 | |
151 emma.delta.REML.LL.w.Z <- function(logdelta, lambda, etas.1, n, t, etas.2.sq ) { | |
152 tq <- length(etas.1) | |
153 nq <- n - t + tq | |
154 delta <- exp(logdelta) | |
155 return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(lambda+delta))+(n-t)*logdelta)) ) | |
156 } | |
157 | |
158 emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) { | |
159 nq <- length(etas) | |
160 delta <- exp(logdelta) | |
161 etasq <- etas*etas | |
162 ldelta <- lambda+delta | |
163 return( 0.5*(nq*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/ldelta)) ) | |
164 } | |
165 | |
166 emma.delta.REML.dLL.w.Z <- function(logdelta, lambda, etas.1, n, t1, etas.2.sq ) { | |
167 t <- t1 | |
168 tq <- length(etas.1) | |
169 nq <- n - t + tq | |
170 delta <- exp(logdelta) | |
171 etasq <- etas.1*etas.1 | |
172 ldelta <- lambda+delta | |
173 return( 0.5*(nq*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/ldelta)+(n-t)/delta)) ) | |
174 } | |
175 | |
176 emma.MLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10, | |
177 esp=1e-10, eig.L = NULL, eig.R = NULL) | |
178 { | |
179 n <- length(y) | |
180 t <- nrow(K) | |
181 q <- ncol(X) | |
182 | |
183 # stopifnot(nrow(K) == t) | |
184 stopifnot(ncol(K) == t) | |
185 stopifnot(nrow(X) == n) | |
186 | |
187 if ( det(crossprod(X,X)) == 0 ) { | |
188 warning("X is singular") | |
189 return (list(ML=0,delta=0,ve=0,vg=0)) | |
190 } | |
191 | |
192 if ( is.null(Z) ) { | |
193 if ( is.null(eig.L) ) { | |
194 eig.L <- emma.eigen.L.wo.Z(K) | |
195 } | |
196 if ( is.null(eig.R) ) { | |
197 eig.R <- emma.eigen.R.wo.Z(K,X) | |
198 } | |
199 etas <- crossprod(eig.R$vectors,y) | |
200 | |
201 | |
202 logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim | |
203 m <- length(logdelta) | |
204 delta <- exp(logdelta) | |
205 Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE) | |
206 Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE) | |
207 Etasq <- matrix(etas*etas,n-q,m) | |
208 LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Xis))) | |
209 dLL <- 0.5*delta*(n*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Xis)) | |
210 | |
211 optlogdelta <- vector(length=0) | |
212 optLL <- vector(length=0) | |
213 if ( dLL[1] < esp ) { | |
214 optlogdelta <- append(optlogdelta, llim) | |
215 optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.R$values,etas,eig.L$values)) | |
216 } | |
217 if ( dLL[m-1] > 0-esp ) { | |
218 optlogdelta <- append(optlogdelta, ulim) | |
219 optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.R$values,etas,eig.L$values)) | |
220 } | |
221 | |
222 for( i in 1:(m-1) ) | |
223 { | |
224 if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
225 { | |
226 r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas, xi=eig.L$values) | |
227 optlogdelta <- append(optlogdelta, r$root) | |
228 optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.R$values, etas, eig.L$values)) | |
229 } | |
230 } | |
231 # optdelta <- exp(optlogdelta) | |
232 } | |
233 else { | |
234 if ( is.null(eig.L) ) { | |
235 eig.L <- emma.eigen.L.w.Z(Z,K) | |
236 } | |
237 if ( is.null(eig.R) ) { | |
238 eig.R <- emma.eigen.R.w.Z(Z,K,X) | |
239 } | |
240 etas <- crossprod(eig.R$vectors,y) | |
241 etas.1 <- etas[1:(t-q)] | |
242 etas.2 <- etas[(t-q+1):(n-q)] | |
243 etas.2.sq <- sum(etas.2*etas.2) | |
244 | |
245 logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim | |
246 | |
247 m <- length(logdelta) | |
248 delta <- exp(logdelta) | |
249 Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE) | |
250 Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE) | |
251 Etasq <- matrix(etas.1*etas.1,t-q,m) | |
252 #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas)) | |
253 dLL <- 0.5*delta*(n*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta)) | |
254 | |
255 optlogdelta <- vector(length=0) | |
256 optLL <- vector(length=0) | |
257 if ( dLL[1] < esp ) { | |
258 optlogdelta <- append(optlogdelta, llim) | |
259 optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq)) | |
260 } | |
261 if ( dLL[m-1] > 0-esp ) { | |
262 optlogdelta <- append(optlogdelta, ulim) | |
263 optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq)) | |
264 } | |
265 | |
266 for( i in 1:(m-1) ) | |
267 { | |
268 if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
269 { | |
270 r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq ) | |
271 optlogdelta <- append(optlogdelta, r$root) | |
272 optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.R$values, etas.1, eig.L$values, n, etas.2.sq )) | |
273 } | |
274 } | |
275 # optdelta <- exp(optlogdelta) | |
276 } | |
277 | |
278 maxdelta <- exp(optlogdelta[which.max(optLL)]) | |
279 maxLL <- max(optLL) | |
280 if ( is.null(Z) ) { | |
281 maxva <- sum(etas*etas/(eig.R$values+maxdelta))/n | |
282 } | |
283 else { | |
284 maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/n | |
285 } | |
286 maxve <- maxva*maxdelta | |
287 | |
288 return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva)) | |
289 } | |
290 | |
291 emma.MLE.noX <- function(y, K, Z=NULL, ngrids=100, llim=-10, ulim=10, | |
292 esp=1e-10, eig.L = NULL) | |
293 { | |
294 n <- length(y) | |
295 t <- nrow(K) | |
296 | |
297 # stopifnot(nrow(K) == t) | |
298 stopifnot(ncol(K) == t) | |
299 | |
300 if ( is.null(Z) ) { | |
301 if ( is.null(eig.L) ) { | |
302 eig.L <- emma.eigen.L.wo.Z(K) | |
303 } | |
304 etas <- crossprod(eig.L$vectors,y) | |
305 | |
306 logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim | |
307 m <- length(logdelta) | |
308 delta <- exp(logdelta) | |
309 Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE) | |
310 Etasq <- matrix(etas*etas,n,m) | |
311 LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Xis)))-colSums(log(Xis))) | |
312 dLL <- 0.5*delta*(n*colSums(Etasq/(Xis*Xis))/colSums(Etasq/Xis)-colSums(1/Xis)) | |
313 | |
314 optlogdelta <- vector(length=0) | |
315 optLL <- vector(length=0) | |
316 #print(dLL) | |
317 if ( dLL[1] < esp ) { | |
318 optlogdelta <- append(optlogdelta, llim) | |
319 optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.L$values,etas,eig.L$values)) | |
320 } | |
321 if ( dLL[m-1] > 0-esp ) { | |
322 optlogdelta <- append(optlogdelta, ulim) | |
323 optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.L$values,etas,eig.L$values)) | |
324 } | |
325 | |
326 for( i in 1:(m-1) ) | |
327 { | |
328 #if ( ( dLL[i]*dLL[i+1] < 0 ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
329 if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
330 { | |
331 r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas=etas, xi=eig.L$values) | |
332 optlogdelta <- append(optlogdelta, r$root) | |
333 optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.L$values, etas, eig.L$values)) | |
334 } | |
335 } | |
336 # optdelta <- exp(optlogdelta) | |
337 } | |
338 else { | |
339 if ( is.null(eig.L) ) { | |
340 eig.L <- emma.eigen.L.w.Z(Z,K) | |
341 } | |
342 etas <- crossprod(eig.L$vectors,y) | |
343 etas.1 <- etas[1:t] | |
344 etas.2 <- etas[(t+1):n] | |
345 etas.2.sq <- sum(etas.2*etas.2) | |
346 | |
347 logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim | |
348 | |
349 m <- length(logdelta) | |
350 delta <- exp(logdelta) | |
351 Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE) | |
352 Etasq <- matrix(etas.1*etas.1,t,m) | |
353 #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas)) | |
354 dLL <- 0.5*delta*(n*(colSums(Etasq/(Xis*Xis))+etas.2.sq/(delta*delta))/(colSums(Etasq/Xis)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta)) | |
355 | |
356 optlogdelta <- vector(length=0) | |
357 optLL <- vector(length=0) | |
358 if ( dLL[1] < esp ) { | |
359 optlogdelta <- append(optlogdelta, llim) | |
360 optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq)) | |
361 } | |
362 if ( dLL[m-1] > 0-esp ) { | |
363 optlogdelta <- append(optlogdelta, ulim) | |
364 optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq)) | |
365 } | |
366 | |
367 for( i in 1:(m-1) ) | |
368 { | |
369 if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
370 { | |
371 r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq ) | |
372 optlogdelta <- append(optlogdelta, r$root) | |
373 optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.L$values, etas.1, eig.L$values, n, etas.2.sq )) | |
374 } | |
375 } | |
376 # optdelta <- exp(optlogdelta) | |
377 } | |
378 | |
379 maxdelta <- exp(optlogdelta[which.max(optLL)]) | |
380 maxLL <- max(optLL) | |
381 if ( is.null(Z) ) { | |
382 maxva <- sum(etas*etas/(eig.L$values+maxdelta))/n | |
383 } | |
384 else { | |
385 maxva <- (sum(etas.1*etas.1/(eig.L$values+maxdelta))+etas.2.sq/maxdelta)/n | |
386 } | |
387 maxve <- maxva*maxdelta | |
388 | |
389 return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva)) | |
390 } | |
391 | |
392 emma.REMLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10, | |
393 esp=1e-10, eig.L = NULL, eig.R = NULL) { | |
394 n <- length(y) | |
395 t <- nrow(K) | |
396 q <- ncol(X) | |
397 | |
398 # stopifnot(nrow(K) == t) | |
399 stopifnot(ncol(K) == t) | |
400 stopifnot(nrow(X) == n) | |
401 | |
402 if ( det(crossprod(X,X)) == 0 ) { | |
403 warning("X is singular") | |
404 return (list(REML=0,delta=0,ve=0,vg=0)) | |
405 } | |
406 | |
407 if ( is.null(Z) ) { | |
408 if ( is.null(eig.R) ) { | |
409 eig.R <- emma.eigen.R.wo.Z(K,X) | |
410 } | |
411 etas <- crossprod(eig.R$vectors,y) | |
412 | |
413 logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim | |
414 m <- length(logdelta) | |
415 delta <- exp(logdelta) | |
416 Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE) | |
417 Etasq <- matrix(etas*etas,n-q,m) | |
418 LL <- 0.5*((n-q)*(log((n-q)/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Lambdas))) | |
419 dLL <- 0.5*delta*((n-q)*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Lambdas)) | |
420 | |
421 optlogdelta <- vector(length=0) | |
422 optLL <- vector(length=0) | |
423 if ( dLL[1] < esp ) { | |
424 optlogdelta <- append(optlogdelta, llim) | |
425 optLL <- append(optLL, emma.delta.REML.LL.wo.Z(llim,eig.R$values,etas)) | |
426 } | |
427 if ( dLL[m-1] > 0-esp ) { | |
428 optlogdelta <- append(optlogdelta, ulim) | |
429 optLL <- append(optLL, emma.delta.REML.LL.wo.Z(ulim,eig.R$values,etas)) | |
430 } | |
431 | |
432 for( i in 1:(m-1) ) | |
433 { | |
434 if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
435 { | |
436 r <- uniroot(emma.delta.REML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas) | |
437 optlogdelta <- append(optlogdelta, r$root) | |
438 optLL <- append(optLL, emma.delta.REML.LL.wo.Z(r$root,eig.R$values, etas)) | |
439 } | |
440 } | |
441 # optdelta <- exp(optlogdelta) | |
442 } | |
443 else { | |
444 if ( is.null(eig.R) ) { | |
445 eig.R <- emma.eigen.R.w.Z(Z,K,X) | |
446 } | |
447 etas <- crossprod(eig.R$vectors,y) | |
448 etas.1 <- etas[1:(t-q)] | |
449 etas.2 <- etas[(t-q+1):(n-q)] | |
450 etas.2.sq <- sum(etas.2*etas.2) | |
451 | |
452 logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim | |
453 m <- length(logdelta) | |
454 delta <- exp(logdelta) | |
455 Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE) | |
456 Etasq <- matrix(etas.1*etas.1,t-q,m) | |
457 dLL <- 0.5*delta*((n-q)*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Lambdas)+(n-t)/delta)) | |
458 | |
459 optlogdelta <- vector(length=0) | |
460 optLL <- vector(length=0) | |
461 if ( dLL[1] < esp ) { | |
462 optlogdelta <- append(optlogdelta, llim) | |
463 optLL <- append(optLL, emma.delta.REML.LL.w.Z(llim,eig.R$values,etas.1,n,t,etas.2.sq)) | |
464 } | |
465 if ( dLL[m-1] > 0-esp ) { | |
466 optlogdelta <- append(optlogdelta, ulim) | |
467 optLL <- append(optLL, emma.delta.REML.LL.w.Z(ulim,eig.R$values,etas.1,n,t,etas.2.sq)) | |
468 } | |
469 | |
470 for( i in 1:(m-1) ) | |
471 { | |
472 if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) | |
473 { | |
474 r <- uniroot(emma.delta.REML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, n=n, t1=t, etas.2.sq = etas.2.sq ) | |
475 optlogdelta <- append(optlogdelta, r$root) | |
476 optLL <- append(optLL, emma.delta.REML.LL.w.Z(r$root,eig.R$values, etas.1, n, t, etas.2.sq )) | |
477 } | |
478 } | |
479 # optdelta <- exp(optlogdelta) | |
480 } | |
481 | |
482 maxdelta <- exp(optlogdelta[which.max(optLL)]) | |
483 maxLL <- max(optLL) | |
484 if ( is.null(Z) ) { | |
485 maxva <- sum(etas*etas/(eig.R$values+maxdelta))/(n-q) | |
486 } | |
487 else { | |
488 maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/(n-q) | |
489 } | |
490 maxve <- maxva*maxdelta | |
491 | |
492 return (list(REML=maxLL,delta=maxdelta,ve=maxve,vg=maxva)) | |
493 } | |
494 | |
495 emma.ML.LRT <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) { | |
496 if ( is.null(dim(ys)) || ncol(ys) == 1 ) { | |
497 ys <- matrix(ys,1,length(ys)) | |
498 } | |
499 if ( is.null(dim(xs)) || ncol(xs) == 1 ) { | |
500 xs <- matrix(xs,1,length(xs)) | |
501 } | |
502 if ( is.null(X0) ) { | |
503 X0 <- matrix(1,ncol(ys),1) | |
504 } | |
505 | |
506 g <- nrow(ys) | |
507 n <- ncol(ys) | |
508 m <- nrow(xs) | |
509 t <- ncol(xs) | |
510 q0 <- ncol(X0) | |
511 q1 <- q0 + 1 | |
512 | |
513 if ( !ponly ) { | |
514 ML1s <- matrix(nrow=m,ncol=g) | |
515 ML0s <- matrix(nrow=m,ncol=g) | |
516 vgs <- matrix(nrow=m,ncol=g) | |
517 ves <- matrix(nrow=m,ncol=g) | |
518 } | |
519 stats <- matrix(nrow=m,ncol=g) | |
520 ps <- matrix(nrow=m,ncol=g) | |
521 ML0 <- vector(length=g) | |
522 | |
523 stopifnot(nrow(K) == t) | |
524 stopifnot(ncol(K) == t) | |
525 stopifnot(nrow(X0) == n) | |
526 | |
527 if ( sum(is.na(ys)) == 0 ) { | |
528 eig.L <- emma.eigen.L(Z,K) | |
529 eig.R0 <- emma.eigen.R(Z,K,X0) | |
530 | |
531 for(i in 1:g) { | |
532 ML0[i] <- emma.MLE(ys[i,],X0,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R0)$ML | |
533 } | |
534 | |
535 x.prev <- vector(length=0) | |
536 | |
537 for(i in 1:m) { | |
538 vids <- !is.na(xs[i,]) | |
539 nv <- sum(vids) | |
540 xv <- xs[i,vids] | |
541 | |
542 if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { | |
543 if (!ponly) { | |
544 stats[i,] <- rep(NA,g) | |
545 vgs[i,] <- rep(NA,g) | |
546 ves[i,] <- rep(NA,g) | |
547 ML1s[i,] <- rep(NA,g) | |
548 ML0s[i,] <- rep(NA,g) | |
549 } | |
550 ps[i,] = rep(1,g) | |
551 } | |
552 else if ( identical(x.prev, xv) ) { | |
553 if ( !ponly ) { | |
554 stats[i,] <- stats[i-1,] | |
555 vgs[i,] <- vgs[i-1,] | |
556 ves[i,] <- ves[i-1,] | |
557 ML1s[i,] <- ML1s[i-1,] | |
558 ML0s[i,] <- ML0s[i-1,] | |
559 } | |
560 ps[i,] <- ps[i-1,] | |
561 } | |
562 else { | |
563 if ( is.null(Z) ) { | |
564 X <- cbind(X0[vids,,drop=FALSE],xs[i,vids]) | |
565 eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X) | |
566 } | |
567 else { | |
568 vrows <- as.logical(rowSums(Z[,vids])) | |
569 nr <- sum(vrows) | |
570 X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids]%*%t(xs[i,vids,drop=FALSE])) | |
571 eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X) | |
572 } | |
573 | |
574 for(j in 1:g) { | |
575 if ( nv == t ) { | |
576 MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1) | |
577 # MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1) | |
578 if (!ponly) { | |
579 ML1s[i,j] <- MLE$ML | |
580 vgs[i,j] <- MLE$vg | |
581 ves[i,j] <- MLE$ve | |
582 } | |
583 stats[i,j] <- 2*(MLE$ML-ML0[j]) | |
584 | |
585 } | |
586 else { | |
587 if ( is.null(Z) ) { | |
588 eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids]) | |
589 MLE0 <- emma.MLE(ys[j,vids],X0[vids,,drop=FALSE],K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0) | |
590 MLE1 <- emma.MLE(ys[j,vids],X,K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0) | |
591 } | |
592 else { | |
593 if ( nr == n ) { | |
594 MLE1 <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L) | |
595 } | |
596 else { | |
597 eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids],K[vids,vids]) | |
598 MLE0 <- emma.MLE(ys[j,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0) | |
599 MLE1 <- emma.MLE(ys[j,vrows],X,K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0) | |
600 } | |
601 } | |
602 if (!ponly) { | |
603 ML1s[i,j] <- MLE1$ML | |
604 ML0s[i,j] <- MLE0$ML | |
605 vgs[i,j] <- MLE1$vg | |
606 ves[i,j] <- MLE1$ve | |
607 } | |
608 stats[i,j] <- 2*(MLE1$ML-MLE0$ML) | |
609 } | |
610 } | |
611 if ( ( nv == t ) && ( !ponly ) ) { | |
612 ML0s[i,] <- ML0 | |
613 } | |
614 ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE) | |
615 } | |
616 } | |
617 } | |
618 else { | |
619 eig.L <- emma.eigen.L(Z,K) | |
620 eig.R0 <- emma.eigen.R(Z,K,X0) | |
621 | |
622 for(i in 1:g) { | |
623 vrows <- !is.na(ys[i,]) | |
624 if ( is.null(Z) ) { | |
625 ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vrows,vrows],NULL,ngrids,llim,ulim,esp)$ML | |
626 } | |
627 else { | |
628 vids <- colSums(Z[vrows,]>0) | |
629 | |
630 ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp)$ML | |
631 } | |
632 } | |
633 | |
634 x.prev <- vector(length=0) | |
635 | |
636 for(i in 1:m) { | |
637 vids <- !is.na(xs[i,]) | |
638 nv <- sum(vids) | |
639 xv <- xs[i,vids] | |
640 | |
641 if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { | |
642 if (!ponly) { | |
643 stats[i,] <- rep(NA,g) | |
644 vgs[i,] <- rep(NA,g) | |
645 ves[i,] <- rep(NA,g) | |
646 ML1s[i,] <- rep(NA,g) | |
647 ML0s[,i] <- rep(NA,g) | |
648 } | |
649 ps[i,] = rep(1,g) | |
650 } | |
651 else if ( identical(x.prev, xv) ) { | |
652 if ( !ponly ) { | |
653 stats[i,] <- stats[i-1,] | |
654 vgs[i,] <- vgs[i-1,] | |
655 ves[i,] <- ves[i-1,] | |
656 ML1s[i,] <- ML1s[i-1,] | |
657 } | |
658 ps[i,] = ps[i-1,] | |
659 } | |
660 else { | |
661 if ( is.null(Z) ) { | |
662 X <- cbind(X0,xs[i,]) | |
663 if ( nv == t ) { | |
664 eig.R1 = emma.eigen.R.wo.Z(K,X) | |
665 } | |
666 } | |
667 else { | |
668 vrows <- as.logical(rowSums(Z[,vids])) | |
669 X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) | |
670 if ( nv == t ) { | |
671 eig.R1 = emma.eigen.R.w.Z(Z,K,X) | |
672 } | |
673 } | |
674 | |
675 for(j in 1:g) { | |
676 # print(j) | |
677 vrows <- !is.na(ys[j,]) | |
678 if ( nv == t ) { | |
679 nr <- sum(vrows) | |
680 if ( is.null(Z) ) { | |
681 if ( nr == n ) { | |
682 MLE <- emma.MLE(ys[j,],X,K,NULL,ngrids,llim,ulim,esp,eig.L,eig.R1) | |
683 } | |
684 else { | |
685 MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vrows,vrows],NULL,ngrids,llim,ulim,esp) | |
686 } | |
687 } | |
688 else { | |
689 if ( nr == n ) { | |
690 MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1) | |
691 } | |
692 else { | |
693 vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) | |
694 MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vtids,vtids],Z[vrows,vtids],ngrids,llim,ulim,esp) | |
695 } | |
696 } | |
697 | |
698 if (!ponly) { | |
699 ML1s[i,j] <- MLE$ML | |
700 vgs[i,j] <- MLE$vg | |
701 ves[i,j] <- MLE$ve | |
702 } | |
703 stats[i,j] <- 2*(MLE$ML-ML0[j]) | |
704 } | |
705 else { | |
706 if ( is.null(Z) ) { | |
707 vtids <- vrows & vids | |
708 eig.L0 <- emma.eigen.L(NULL,K[vtids,vtids]) | |
709 MLE0 <- emma.MLE(ys[j,vtids],X0[vtids,,drop=FALSE],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0) | |
710 MLE1 <- emma.MLE(ys[j,vtids],X[vtids,],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0) | |
711 } | |
712 else { | |
713 vtids <- as.logical(colSums(Z[vrows,])) & vids | |
714 vtrows <- vrows & as.logical(rowSums(Z[,vids])) | |
715 eig.L0 <- emma.eigen.L(Z[vtrows,vtids],K[vtids,vtids]) | |
716 MLE0 <- emma.MLE(ys[j,vtrows],X0[vtrows,,drop=FALSE],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0) | |
717 MLE1 <- emma.MLE(ys[j,vtrows],X[vtrows,],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0) | |
718 } | |
719 if (!ponly) { | |
720 ML1s[i,j] <- MLE1$ML | |
721 vgs[i,j] <- MLE1$vg | |
722 ves[i,j] <- MLE1$ve | |
723 ML0s[i,j] <- MLE0$ML | |
724 } | |
725 stats[i,j] <- 2*(MLE1$ML-MLE0$ML) | |
726 } | |
727 } | |
728 if ( ( nv == t ) && ( !ponly ) ) { | |
729 ML0s[i,] <- ML0 | |
730 } | |
731 ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE) | |
732 } | |
733 } | |
734 } | |
735 if ( ponly ) { | |
736 return (ps) | |
737 } | |
738 else { | |
739 return (list(ps=ps,ML1s=ML1s,ML0s=ML0s,stats=stats,vgs=vgs,ves=ves)) | |
740 } | |
741 } | |
742 | |
743 emma.test <- function(ys, xs, K, Z=NULL, x0s = NULL, X0 = NULL, dfxs = 1, dfx0s = 1, use.MLE = FALSE, use.LRT = FALSE, ngrids = 100, llim = -10, ulim = 10, esp=1e-10, ponly = FALSE) | |
744 { | |
745 stopifnot (dfxs > 0) | |
746 | |
747 if ( is.null(dim(ys)) || ncol(ys) == 1 ) { | |
748 ys <- matrix(ys,1,length(ys)) | |
749 } | |
750 | |
751 if ( is.null(dim(xs)) || ncol(xs) == 1 ) { | |
752 xs <- matrix(xs,1,length(xs)) | |
753 } | |
754 nx <- nrow(xs)/dfxs | |
755 | |
756 if ( is.null(x0s) ) { | |
757 dfx0s = 0 | |
758 x0s <- matrix(NA,0,ncol(xs)) | |
759 } | |
760 # X0 automatically contains intercept. If no intercept is to be used, | |
761 # X0 should be matrix(nrow=ncol(ys),ncol=0) | |
762 if ( is.null(X0) ) { | |
763 X0 <- matrix(1,ncol(ys),1) | |
764 } | |
765 | |
766 stopifnot(Z == NULL) # The case where Z is not null is not implemented | |
767 | |
768 ny <- nrow(ys) | |
769 iy <- ncol(ys) | |
770 ix <- ncol(xs) | |
771 | |
772 stopifnot(nrow(K) == ix) | |
773 stopifnot(ncol(K) == ix) | |
774 stopifnot(nrow(X0) == iy) | |
775 | |
776 if ( !ponly ) { | |
777 LLs <- matrix(nrow=m,ncol=g) | |
778 vgs <- matrix(nrow=m,ncol=g) | |
779 ves <- matrix(nrow=m,ncol=g) | |
780 } | |
781 dfs <- matrix(nrow=m,ncol=g) | |
782 stats <- matrix(nrow=m,ncol=g) | |
783 ps <- matrix(nrow=m,ncol=g) | |
784 | |
785 # The case with no missing phenotypes | |
786 if ( sum(is.na(ys)) == 0 ) { | |
787 if ( ( use.MLE ) || ( !use.LRT ) ) { | |
788 eig.L0 <- emma.eigen.L(Z,K) | |
789 } | |
790 if ( dfx0s == 0 ) { | |
791 eig.R0 <- emma.eigen.R(Z,K,X0) | |
792 } | |
793 x.prev <- NULL | |
794 | |
795 for(i in 1:ix) { | |
796 x1 <- t(xs[(dfxs*(i-1)+1):(dfxs*i),,drop=FALSE]) | |
797 if ( dfxs0 == 0 ) { | |
798 x0 <- X0 | |
799 } | |
800 else { | |
801 x0 <- cbind(t(x0s[(dfx0s*(i-1)+1):(dfx0s*i),,drop=FALSE]),X0) | |
802 } | |
803 x <- cbind(x1,x0) | |
804 xvids <- rowSums(is.na(x) == 0) | |
805 nxv <- sum(xvids) | |
806 xv <- x[xvids,,drop=FALSE] | |
807 Kv <- K[xvids,xvids,drop=FALSE] | |
808 yv <- ys[j,xvids] | |
809 | |
810 if ( identical(x.prev, xv) ) { | |
811 if ( !ponly ) { | |
812 vgs[i,] <- vgs[i-1,] | |
813 ves[i,] <- ves[i-1,] | |
814 dfs[i,] <- dfs[i-1,] | |
815 REMLs[i,] <- REMLs[i-1,] | |
816 stats[i,] <- stats[i-1,] | |
817 } | |
818 ps[i,] <- ps[i-1,] | |
819 } | |
820 else { | |
821 eig.R1 = emma.eigen.R.wo.Z(Kv,xv) | |
822 | |
823 for(j in 1:iy) { | |
824 if ( ( use.MLE ) || ( !use.LRT ) ) { | |
825 if ( nxv < t ) { | |
826 # NOTE: this complexity can be improved by avoiding eigen computation for identical missing patterns | |
827 eig.L0v <- emma.eigen.L.wo.Z(Kv) | |
828 } | |
829 else { | |
830 eig.L0v <- eig.L0 | |
831 } | |
832 } | |
833 | |
834 if ( use.MLE ) { | |
835 MLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1) | |
836 stop("Not implemented yet") | |
837 } | |
838 else { | |
839 REMLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1) | |
840 if ( use.LRT ) { | |
841 stop("Not implemented yet") | |
842 } | |
843 else { | |
844 U <- eig.L0v$vectors * matrix(sqrt(1/(eig.L0v$values+REMLE$delta)),t,t,byrow=TRUE) | |
845 dfs[i,j] <- length(eig.R1$values) | |
846 yt <- crossprod(U,yv) | |
847 xt <- crossprod(U,xv) | |
848 ixx <- solve(crossprod(xt,xt)) | |
849 beta <- ixx%*%crossprod(xt,yt) | |
850 if ( dfxs == 1 ) { | |
851 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
852 } | |
853 else { | |
854 model.m <- c(rep(1,dfxs),rep(0,ncol(xv)-dfxs)) | |
855 stats[i,j] <- | |
856 crossprod(crossprod(solve(crossprod(crossprod(iXX,model.m), | |
857 model.m)), | |
858 model.m*beta),model.m*beta) | |
859 | |
860 } | |
861 if ( !ponly ) { | |
862 vgs[i,j] <- REMLE$vg | |
863 ves[i,j] <- REMLE$ve | |
864 REMLs[i,j] <- REMLE$REML | |
865 } | |
866 } | |
867 } | |
868 } | |
869 if ( dfxs == 1 ) { | |
870 ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) | |
871 } | |
872 else { | |
873 ps[i,] <- pf(abs(stats[i,]),dfs[i,],lower.tail=FALSE) | |
874 } | |
875 } | |
876 } | |
877 } | |
878 # The case with missing genotypes - not implemented yet | |
879 else { | |
880 stop("Not implemented yet") | |
881 eig.L <- emma.eigen.L(Z,K) | |
882 eig.R0 <- emma.eigen.R(Z,K,X0) | |
883 | |
884 x.prev <- vector(length=0) | |
885 | |
886 for(i in 1:m) { | |
887 vids <- !is.na(xs[i,]) | |
888 nv <- sum(vids) | |
889 xv <- xs[i,vids] | |
890 | |
891 if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { | |
892 if (!ponly) { | |
893 vgs[i,] <- rep(NA,g) | |
894 ves[i,] <- rep(NA,g) | |
895 REMLs[i,] <- rep(NA,g) | |
896 dfs[i,] <- rep(NA,g) | |
897 } | |
898 ps[i,] = rep(1,g) | |
899 } | |
900 else if ( identical(x.prev, xv) ) { | |
901 if ( !ponly ) { | |
902 stats[i,] <- stats[i-1,] | |
903 vgs[i,] <- vgs[i-1,] | |
904 ves[i,] <- ves[i-1,] | |
905 REMLs[i,] <- REMLs[i-1,] | |
906 dfs[i,] <- dfs[i-1,] | |
907 } | |
908 ps[i,] = ps[i-1,] | |
909 } | |
910 else { | |
911 if ( is.null(Z) ) { | |
912 X <- cbind(X0,xs[i,]) | |
913 if ( nv == t ) { | |
914 eig.R1 = emma.eigen.R.wo.Z(K,X) | |
915 } | |
916 } | |
917 else { | |
918 vrows <- as.logical(rowSums(Z[,vids,drop=FALSE])) | |
919 X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) | |
920 if ( nv == t ) { | |
921 eig.R1 = emma.eigen.R.w.Z(Z,K,X) | |
922 } | |
923 } | |
924 | |
925 for(j in 1:g) { | |
926 vrows <- !is.na(ys[j,]) | |
927 if ( nv == t ) { | |
928 yv <- ys[j,vrows] | |
929 nr <- sum(vrows) | |
930 if ( is.null(Z) ) { | |
931 if ( nr == n ) { | |
932 REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1) | |
933 U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) | |
934 } | |
935 else { | |
936 eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE]) | |
937 REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp) | |
938 U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) | |
939 } | |
940 dfs[i,j] <- nr-q1 | |
941 } | |
942 else { | |
943 if ( nr == n ) { | |
944 REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1) | |
945 U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) | |
946 } | |
947 else { | |
948 vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) | |
949 eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) | |
950 REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) | |
951 U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) | |
952 } | |
953 dfs[i,j] <- nr-q1 | |
954 } | |
955 | |
956 yt <- crossprod(U,yv) | |
957 Xt <- crossprod(U,X[vrows,,drop=FALSE]) | |
958 iXX <- solve(crossprod(Xt,Xt)) | |
959 beta <- iXX%*%crossprod(Xt,yt) | |
960 if ( !ponly ) { | |
961 vgs[i,j] <- REMLE$vg | |
962 ves[i,j] <- REMLE$ve | |
963 REMLs[i,j] <- REMLE$REML | |
964 } | |
965 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
966 } | |
967 else { | |
968 if ( is.null(Z) ) { | |
969 vtids <- vrows & vids | |
970 eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE]) | |
971 yv <- ys[j,vtids] | |
972 nr <- sum(vtids) | |
973 REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp) | |
974 U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) | |
975 Xt <- crossprod(U,X[vtids,,drop=FALSE]) | |
976 dfs[i,j] <- nr-q1 | |
977 } | |
978 else { | |
979 vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids | |
980 vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE])) | |
981 eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) | |
982 yv <- ys[j,vtrows] | |
983 nr <- sum(vtrows) | |
984 REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) | |
985 U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) | |
986 Xt <- crossprod(U,X[vtrows,,drop=FALSE]) | |
987 dfs[i,j] <- nr-q1 | |
988 } | |
989 yt <- crossprod(U,yv) | |
990 iXX <- solve(crossprod(Xt,Xt)) | |
991 beta <- iXX%*%crossprod(Xt,yt) | |
992 if ( !ponly ) { | |
993 vgs[i,j] <- REMLE$vg | |
994 ves[i,j] <- REMLE$ve | |
995 REMLs[i,j] <- REMLE$REML | |
996 } | |
997 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
998 | |
999 } | |
1000 } | |
1001 ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) | |
1002 } | |
1003 } | |
1004 } | |
1005 if ( ponly ) { | |
1006 return (ps) | |
1007 } | |
1008 else { | |
1009 return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves)) | |
1010 } | |
1011 } | |
1012 | |
1013 emma.REML.t <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) { | |
1014 if ( is.null(dim(ys)) || ncol(ys) == 1 ) { | |
1015 ys <- matrix(ys,1,length(ys)) | |
1016 } | |
1017 if ( is.null(dim(xs)) || ncol(xs) == 1 ) { | |
1018 xs <- matrix(xs,1,length(xs)) | |
1019 } | |
1020 if ( is.null(X0) ) { | |
1021 X0 <- matrix(1,ncol(ys),1) | |
1022 } | |
1023 | |
1024 g <- nrow(ys) | |
1025 n <- ncol(ys) | |
1026 m <- nrow(xs) | |
1027 t <- ncol(xs) | |
1028 q0 <- ncol(X0) | |
1029 q1 <- q0 + 1 | |
1030 | |
1031 stopifnot(nrow(K) == t) | |
1032 stopifnot(ncol(K) == t) | |
1033 stopifnot(nrow(X0) == n) | |
1034 | |
1035 if ( !ponly ) { | |
1036 REMLs <- matrix(nrow=m,ncol=g) | |
1037 vgs <- matrix(nrow=m,ncol=g) | |
1038 ves <- matrix(nrow=m,ncol=g) | |
1039 } | |
1040 dfs <- matrix(nrow=m,ncol=g) | |
1041 stats <- matrix(nrow=m,ncol=g) | |
1042 ps <- matrix(nrow=m,ncol=g) | |
1043 | |
1044 if ( sum(is.na(ys)) == 0 ) { | |
1045 eig.L <- emma.eigen.L(Z,K) | |
1046 | |
1047 x.prev <- vector(length=0) | |
1048 | |
1049 for(i in 1:m) { | |
1050 vids <- !is.na(xs[i,]) | |
1051 nv <- sum(vids) | |
1052 xv <- xs[i,vids] | |
1053 | |
1054 if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { | |
1055 if ( !ponly ) { | |
1056 vgs[i,] <- rep(NA,g) | |
1057 ves[i,] <- rep(NA,g) | |
1058 dfs[i,] <- rep(NA,g) | |
1059 REMLs[i,] <- rep(NA,g) | |
1060 stats[i,] <- rep(NA,g) | |
1061 } | |
1062 ps[i,] = rep(1,g) | |
1063 | |
1064 } | |
1065 else if ( identical(x.prev, xv) ) { | |
1066 if ( !ponly ) { | |
1067 vgs[i,] <- vgs[i-1,] | |
1068 ves[i,] <- ves[i-1,] | |
1069 dfs[i,] <- dfs[i-1,] | |
1070 REMLs[i,] <- REMLs[i-1,] | |
1071 stats[i,] <- stats[i-1,] | |
1072 } | |
1073 ps[i,] <- ps[i-1,] | |
1074 } | |
1075 else { | |
1076 if ( is.null(Z) ) { | |
1077 X <- cbind(X0[vids,,drop=FALSE],xs[i,vids]) | |
1078 eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X) | |
1079 } | |
1080 else { | |
1081 vrows <- as.logical(rowSums(Z[,vids])) | |
1082 X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) | |
1083 eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X) | |
1084 } | |
1085 | |
1086 for(j in 1:g) { | |
1087 if ( nv == t ) { | |
1088 REMLE <- emma.REMLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.R1) | |
1089 if ( is.null(Z) ) { | |
1090 U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),t,t,byrow=TRUE) | |
1091 dfs[i,j] <- nv - q1 | |
1092 } | |
1093 else { | |
1094 U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) | |
1095 dfs[i,j] <- n - q1 | |
1096 } | |
1097 yt <- crossprod(U,ys[j,]) | |
1098 Xt <- crossprod(U,X) | |
1099 iXX <- solve(crossprod(Xt,Xt)) | |
1100 beta <- iXX%*%crossprod(Xt,yt) | |
1101 | |
1102 if ( !ponly ) { | |
1103 vgs[i,j] <- REMLE$vg | |
1104 ves[i,j] <- REMLE$ve | |
1105 REMLs[i,j] <- REMLE$REML | |
1106 } | |
1107 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
1108 } | |
1109 else { | |
1110 if ( is.null(Z) ) { | |
1111 eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids]) | |
1112 nr <- sum(vids) | |
1113 yv <- ys[j,vids] | |
1114 REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],NULL,ngrids,llim,ulim,esp,eig.R1) | |
1115 U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) | |
1116 dfs[i,j] <- nr - q1 | |
1117 } | |
1118 else { | |
1119 eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids,drop=FALSE],K[vids,vids]) | |
1120 yv <- ys[j,vrows] | |
1121 nr <- sum(vrows) | |
1122 tv <- sum(vids) | |
1123 REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],Z[vrows,vids,drop=FALSE],ngrids,llim,ulim,esp,eig.R1) | |
1124 U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-tv)),nr,nr,byrow=TRUE) | |
1125 dfs[i,j] <- nr - q1 | |
1126 } | |
1127 yt <- crossprod(U,yv) | |
1128 Xt <- crossprod(U,X) | |
1129 iXX <- solve(crossprod(Xt,Xt)) | |
1130 beta <- iXX%*%crossprod(Xt,yt) | |
1131 if (!ponly) { | |
1132 vgs[i,j] <- REMLE$vg | |
1133 ves[i,j] <- REMLE$ve | |
1134 REMLs[i,j] <- REMLE$REML | |
1135 } | |
1136 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
1137 } | |
1138 } | |
1139 ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) | |
1140 } | |
1141 } | |
1142 } | |
1143 else { | |
1144 eig.L <- emma.eigen.L(Z,K) | |
1145 eig.R0 <- emma.eigen.R(Z,K,X0) | |
1146 | |
1147 x.prev <- vector(length=0) | |
1148 | |
1149 for(i in 1:m) { | |
1150 vids <- !is.na(xs[i,]) | |
1151 nv <- sum(vids) | |
1152 xv <- xs[i,vids] | |
1153 | |
1154 if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) { | |
1155 if (!ponly) { | |
1156 vgs[i,] <- rep(NA,g) | |
1157 ves[i,] <- rep(NA,g) | |
1158 REMLs[i,] <- rep(NA,g) | |
1159 dfs[i,] <- rep(NA,g) | |
1160 } | |
1161 ps[i,] = rep(1,g) | |
1162 } | |
1163 else if ( identical(x.prev, xv) ) { | |
1164 if ( !ponly ) { | |
1165 stats[i,] <- stats[i-1,] | |
1166 vgs[i,] <- vgs[i-1,] | |
1167 ves[i,] <- ves[i-1,] | |
1168 REMLs[i,] <- REMLs[i-1,] | |
1169 dfs[i,] <- dfs[i-1,] | |
1170 } | |
1171 ps[i,] = ps[i-1,] | |
1172 } | |
1173 else { | |
1174 if ( is.null(Z) ) { | |
1175 X <- cbind(X0,xs[i,]) | |
1176 if ( nv == t ) { | |
1177 eig.R1 = emma.eigen.R.wo.Z(K,X) | |
1178 } | |
1179 } | |
1180 else { | |
1181 vrows <- as.logical(rowSums(Z[,vids,drop=FALSE])) | |
1182 X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE])) | |
1183 if ( nv == t ) { | |
1184 eig.R1 = emma.eigen.R.w.Z(Z,K,X) | |
1185 } | |
1186 } | |
1187 | |
1188 for(j in 1:g) { | |
1189 vrows <- !is.na(ys[j,]) | |
1190 if ( nv == t ) { | |
1191 yv <- ys[j,vrows] | |
1192 nr <- sum(vrows) | |
1193 if ( is.null(Z) ) { | |
1194 if ( nr == n ) { | |
1195 REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1) | |
1196 U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) | |
1197 } | |
1198 else { | |
1199 eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE]) | |
1200 REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp) | |
1201 U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) | |
1202 } | |
1203 dfs[i,j] <- nr-q1 | |
1204 } | |
1205 else { | |
1206 if ( nr == n ) { | |
1207 REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1) | |
1208 U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) | |
1209 } | |
1210 else { | |
1211 vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) | |
1212 eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) | |
1213 REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) | |
1214 U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) | |
1215 } | |
1216 dfs[i,j] <- nr-q1 | |
1217 } | |
1218 | |
1219 yt <- crossprod(U,yv) | |
1220 Xt <- crossprod(U,X[vrows,,drop=FALSE]) | |
1221 iXX <- solve(crossprod(Xt,Xt)) | |
1222 beta <- iXX%*%crossprod(Xt,yt) | |
1223 if ( !ponly ) { | |
1224 vgs[i,j] <- REMLE$vg | |
1225 ves[i,j] <- REMLE$ve | |
1226 REMLs[i,j] <- REMLE$REML | |
1227 } | |
1228 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
1229 } | |
1230 else { | |
1231 if ( is.null(Z) ) { | |
1232 vtids <- vrows & vids | |
1233 eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE]) | |
1234 yv <- ys[j,vtids] | |
1235 nr <- sum(vtids) | |
1236 REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp) | |
1237 U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE) | |
1238 Xt <- crossprod(U,X[vtids,,drop=FALSE]) | |
1239 dfs[i,j] <- nr-q1 | |
1240 } | |
1241 else { | |
1242 vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids | |
1243 vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE])) | |
1244 eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE]) | |
1245 yv <- ys[j,vtrows] | |
1246 nr <- sum(vtrows) | |
1247 REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp) | |
1248 U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE) | |
1249 Xt <- crossprod(U,X[vtrows,,drop=FALSE]) | |
1250 dfs[i,j] <- nr-q1 | |
1251 } | |
1252 yt <- crossprod(U,yv) | |
1253 iXX <- solve(crossprod(Xt,Xt)) | |
1254 beta <- iXX%*%crossprod(Xt,yt) | |
1255 if ( !ponly ) { | |
1256 vgs[i,j] <- REMLE$vg | |
1257 ves[i,j] <- REMLE$ve | |
1258 REMLs[i,j] <- REMLE$REML | |
1259 } | |
1260 stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg) | |
1261 | |
1262 } | |
1263 } | |
1264 ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) | |
1265 } | |
1266 } | |
1267 } | |
1268 if ( ponly ) { | |
1269 return (ps) | |
1270 } | |
1271 else { | |
1272 return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves)) | |
1273 } | |
1274 } |