swh:1:snp:3a44eb759780145deea094ac2a25c5049546a085
Tip revision: 65415c7011cbb1bdd6d5b9de5895d2061d3fb471 authored by Han Lin Shang on 01 June 2023, 16:00:02 UTC
version 6.2
version 6.2
Tip revision: 65415c7
MFPCA.R
MFPCA <-
function(y, M=NULL, J=NULL, N=NULL)
{
### Calculate overall mean function and the deviation
### from the overall mean to country specific mean functions
#### Estimate mu ###
mu<- colMeans(y)
### Estimate eta ###
eta <- matrix(0, M, N)
for(m in 1:M) {
eta[ m, ] <- colMeans(y[ (1:J)+(m-1)*J, ]) - mu
}
resd <- matrix(0, nrow=M*J, ncol=N)
for(m in 1:M) {
resd[ (1:J)+(m-1)*J, ] <- t ( t( y[ (1:J)+(m-1)*J, ] ) - (mu + eta[m,]) )
}
### Estimate Rj
Rj<-matrix(0, J, N)
for (j in 1:J){
Rj[j,]<-colMeans(resd[(0:(M-1)*J) + j,])
}
### Estimate Uij
R<-matrix( rep( t( Rj ) , M ) , ncol = ncol(Rj) , byrow = TRUE )
Uij<-resd - R
### Estimate the covariance functions for time trend and functional pattern
G1 <- cov(eta)
G2 <- long_run_covariance_estimation(t(Rj))
G3 <- cov(Uij)
### Estimate eigen values and eigen functions at two levels by calling the
### eigen function (in R "base" package) on discretized covariance matrices.
e1 <- eigen(G1)
e2 <- eigen(G2)
e3 <- eigen(G3)
### transform the eigenvalues of the covariance matrices into eigenvalues
### of original functions
fpca1.value <- e1$values
fpca2.value <- e2$values
fpca3.value <- e3$values
### Keep only non-negative eigenvalues
fpca1.value <- ifelse(fpca1.value>=0, fpca1.value, 0)
fpca2.value <- ifelse(fpca2.value>=0, fpca2.value, 0)
fpca3.value <- ifelse(fpca3.value>=0, fpca3.value, 0)
### Calculate the percentage of variance that are explained by the components
percent1 <- (fpca1.value)/sum(fpca1.value)
percent2 <- (fpca2.value)/sum(fpca2.value)
percent3 <- (fpca3.value)/sum(fpca3.value)
### Calculate the ratio of first eigenvalue and the rest eigenvalues
ratio1<- fpca1.value[1]/fpca1.value
ratio2<- fpca2.value[1]/fpca2.value
ratio3<- fpca3.value[1]/fpca3.value
### Decide the number of components that are kept at level 1 and 2. The general
### rule is to stop at the component where the cumulative percentage of variance
### explained is greater than 90%.
K1 <- max(min(which(cumsum(percent1) > 0.9)), min(which(ratio1>sqrt(M)/log10(M)))-1, 2)
K2 <- max(min(which(cumsum(percent2) > 0.9)), min(which(ratio2>sqrt(J)/log10(J)))-1, 2)
K3 <- max(min(which(cumsum(percent3) > 0.9)), min(which(ratio3>sqrt(M*J)/log10(M*J)))-1, 2)
### estimate eigen vectors for discretized covariance matrices and
### transform them into norm one eigenfunctions
fpca1.vectors <- e1$vectors[, 1:K1]
fpca2.vectors <- e2$vectors[, 1:K2]
fpca3.vectors <- e3$vectors[, 1:K3]
for(i in 1:K1) {
v2 <- fpca1.vectors[,i]
tempsign <- sum(v2)
fpca1.vectors[,i] <- ifelse(tempsign<0, -1,1) * v2
}
for(i in 1:K2) {
v2 <- fpca2.vectors[,i]
tempsign <- sum(v2)
fpca2.vectors[,i] <- ifelse(tempsign<0, -1,1) * v2
}
for(i in 1:K3) {
v2 <- fpca3.vectors[,i]
tempsign <- sum(v2)
fpca3.vectors[,i] <- ifelse(tempsign<0, -1,1) * v2
}
s1 <- eta%*%fpca1.vectors
s2 <- Rj%*%fpca2.vectors
s3 <- Uij%*%fpca3.vectors
### Return the results from multilevel FPCA as a list
return( list(K1=K1, K2=K2,K3=K3, lambda1=fpca1.value, lambda2=fpca2.value, lambda3 = fpca3.value, phi1 = fpca1.vectors, phi2=fpca2.vectors,
phi3=fpca3.vectors, scores1=s1, scores2=s2, scores3 = s3, mu=mu, eta=t(eta), Rj =Rj, Uij = Uij) )
}