Mercurial > repos > dereeper > mlmm

comparison mlmm/source_library/emma.mlmm.r @ 0:6b7107812931 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Thu, 02 Jul 2015 05:42:38 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6b7107812931
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 }