Mercurial > repos > dereeper > mlmm

comparison source_library/emma.r @ 1:380b364980f9 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 }