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