|
0
|
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 }
|
