[程序汇编]Markov Models using R and Matlab

1. Examples
   
2. # Set the number of threads
   
3. nth <- setThreadsTP(2)
   
4. # Create survTP object
   
5. data(heartTP)
   
6. heartTP_obj <- with(heartTP, survTP(time1, event1, Stime, event))
   
7. 48 transPAJ
   
8. # Compute transition probabilities
   
9. transPAJ(object=heartTP_obj, s=33, t=412)
   
10. # Compute transition probabilities with confidence band
    
11. transPAJ(object=heartTP_obj, s=33, t=412, conf=TRUE, conf.level=0.9,
    
12. method.boot="percentile")
    
13. # Restore the number of threads
    
14. setThreadsTP(nth)

复制代码

https://cran.r-project.org/web/packages/TPmsm/TPmsm.pdf

1. Examples
   
2. # draw a sample for 1000 units and only one response variable and 5 time occasions
   
3. k1 = 2; k2 = 3
   
4. la = rep(1/k1,k1)
   
5. Piv = matrix(1/k2,k2,k1)
   
6. Pi = array(0,c(k2,k2,k1))
   
7. Pi[,,1] = diag(k2)
   
8. Pi[,,2] = 1/k2
   
9. Psi = cbind(c(0.6,0.3,0.1),c(0.1,0.3,0.6),c(0.3,0.6,0.1))
   
10. out = draw_lm_mixed(la,Piv,Pi,Psi,n=1000,TT=5)

复制代码

https://cran.r-project.org/web/packages/LMest/LMest.pdf

Estimate Basic LM Model

1. Examples
   
2. # Example of drug consumption data
   
3. # load data
   
4. data(data_drug)
   
5. data_drug = as.matrix(data_drug)
   
6. S = data_drug[,1:5]-1
   
7. yv = data_drug[,6]
   
8. # fit of the Basic LC model
   
9. k = 3
   
10. out = est_lm_basic(S,yv,k,mod=1)
    
11. ## Not run:
    
12. # Example based on criminal data
    
13. # load criminal data
    
14. data(data_criminal_sim)
    
15. est_lm_cov_latent 17
    
16. out = long2wide(data_criminal_sim,"id","time","sex",
    
17. c("y1","y2","y3","y4","y5","y6","y7","y8","y9","y10"),aggr=T,full=999)
    
18. XX = out$XX
    
19. YY = out$YY
    
20. freq = out$freq
    
21. # fit basic LM model with increasing number of states to select the most suitable
    
22. Res0 = vector("list",7)
    
23. for(k in 1:7){
    
24. Res0[[k]] = est_lm_basic(YY,freq,k,mod=1,tol=10^-4)
    
25. save(list = ls(),file="example_criminal_temp.RData")
    
26. }
    
27. out1 = Res0[[6]]
    
28. ## End(Not run)

复制代码

https://cran.r-project.org/web/packages/LMest/LMest.pdf

Estimate LM Model with Covariates in the Latent Model

1. Examples
   
2. ## Not run:
   
3. # Example based on self-rated health status (SRHS) data
   
4. # load SRHS data
   
5. data(data_SRHS_long)
   
6. dataSRHS = data_SRHS_long
   
7. TT = 8
   
8. head(dataSRHS)
   
9. res = long2matrices(dataSRHS$id,X=cbind(dataSRHS$gender-1,
   
10. dataSRHS$race==2|dataSRHS$race==3, dataSRHS$education==4,
    
11. dataSRHS$education==5,dataSRHS$age-50,(dataSRHS$age-50)^2/100),
    
12. Y=dataSRHS$srhs)
    
13. # matrix of responses (with ordered categories from 0 to 4)
    
14. S = 5-res$YY
    
15. n = dim(S)[1]
    
16. # matrix of covariates (for the first and the following occasions)
    
17. # colums are: gender,race,educational level (2 columns),age,age^2)
    
18. X1 =res$XX[,1,]
    
19. X2 =res$XX[,2:TT,]
    
20. # estimate the model
    
21. est2f = est_lm_cov_latent(S,X1,X2,k=2,output=TRUE,out_se=TRUE)
    
22. # average transition probability matrix
    
23. PI = round(apply(est2f$PI[,,,2:TT],c(1,2),mean),4)
    
24. # Transition probability matrix for white females with high educational level
    
25. ind1 = (X1[,1]==1 & X1[,2]==0 & X1[,4]==1)
    
26. PI1 = round(apply(est2f$PI[,,ind1,2:TT],c(1,2),mean),4)
    
27. # Transition probability matrix for non-white male, low educational level
    
28. ind2 = (X1[,1]==0 & X1[,2]==1& X1[,3]==0 & X1[,4]==0)
    
29. PI2 = round(apply(est2f$PI[,,ind2,2:TT],c(1,2),mean),4)
    
30. ## End(Not run)

复制代码

https://cran.r-project.org/web/packages/LMest/LMest.pdf

Estimate LM Model with Covariates in the Measurement Model

1. Examples
   
2. ## Not run:
   
3. # Example based on self-rated health status (SRHS) data
   
4. # load SRHS data
   
5. data(data_SRHS_long)
   
6. dataSRHS = data_SRHS_long
   
7. head(dataSRHS)
   
8. res = long2matrices(dataSRHS$id,X=cbind(dataSRHS$gender-1,
   
9. dataSRHS$race==2|dataSRHS$race==3,dataSRHS$education==4,
   
10. dataSRHS$education==5,dataSRHS$age-50,(dataSRHS$age-50)^2/100),
    
11. Y=dataSRHS$srhs)
    
12. X =res$XX
    
13. S = 5-res$YY
    
14. # *** fit stationary LM model
    
15. res0 = vector("list",10); tol = 10^-6;
    
16. for(k in 1:10){
    
17. res0[[k]] = est_lm_cov_manifest(S,X,k,1,mod=0,tol)
    
18. save.image("example_SRHS.RData")
    
19. }
    
20. # *** fit the mixture latent auto-regressive model
    
21. tol = 0.005
    
22. res = vector("list",4)
    
23. k=1
    
24. q = 51
    
25. res[[k]]=est_lm_cov_manifest(S,X,k,q,mod=1,tol,output=T)
    
26. for(k in 2:4) res[[k]]=est_lm_cov_manifest(S,X,k,q=61,mod=1,tol,output=T)
    
27. ## End(Not run)

复制代码

https://cran.r-project.org/web/packages/LMest/LMest.pdf

Markov Chain Monte Carlo for the Hidden Markov Random-effects Model

1. Examples
   
2. ## Not run:
   
3. ## data generating
   
4. set.seed(1977)
   
5. Q <- 3
   
6. true.beta1 <- c(1, 1, 1) ; true.beta2 <- c(-1, -1, -1)
   
7. true.sigma2 <- c(2, 5); true.D1 <- diag(.5, Q); true.D2 <- diag(2.5, Q)
   
8. N=30; T=100;
   
9. NT <- N*T
   
10. x1 <- runif(NT, 1, 2)
    
11. x2 <- runif(NT, 1, 2)
    
12. X <- cbind(1, x1, x2); W <- X; y <- rep(NA, NT)
    
13. ## true break numbers are one and at the center
    
14. break.point = rep(T/2, N); break.sigma=c(rep(1, N));
    
15. break.list <- rep(1, N)
    
16. id <- rep(1:N, each=NT/N)
    
17. K <- ncol(X);
    
18. ruler <- c(1:T)
    
19. ## compute the weight for the break
    
20. W.mat <- matrix(NA, T, N)
    
21. for (i in 1:N){
    
22. W.mat[, i] <- pnorm((ruler-break.point[i])/break.sigma[i])
    
23. }
    
24. Weight <- as.vector(W.mat)
    
25. ## data generating by weighting two means and variances
    
26. j = 1
    
27. for (i in 1:N){
    
28. Xi <- X[j:(j+T-1), ]
    
29. Wi <- W[j:(j+T-1), ]
    
30. true.V1 <- true.sigma2[1]*diag(T) + Wi%*%true.D1%*%t(Wi)
    
31. true.V2 <- true.sigma2[2]*diag(T) + Wi%*%true.D2%*%t(Wi)
    
32. true.mean1 <- Xi%*%true.beta1
    
33. true.mean2 <- Xi%*%true.beta2
    
34. weight <- Weight[j:(j+T-1)]
    
35. y[j:(j+T-1)] <- (1-weight)*true.mean1 + (1-weight)*chol(true.V1)%*%rnorm(T) +
    
36. weight*true.mean2 + weight*chol(true.V2)%*%rnorm(T)
    
37. j <- j + T
    
38. }
    
39. ## model fitting
    
40. subject.id <- c(rep(1:N, each=T))
    
41. time.id <- c(rep(1:T, N))
    
42. ## model fitting
    
43. G <- 100
    
44. b0 <- rep(0, K) ; B0 <- solve(diag(100, K))
    
45. c0 <- 2; d0 <- 2
    
46. r0 <- 5; R0 <- diag(c(1, 0.1, 0.1))
    
47. subject.id <- c(rep(1:N, each=T))
    
48. time.id <- c(rep(1:T, N))
    
49. out1 <- HMMpanelRE(subject.id, time.id, y, X, W, m=1,
    
50. mcmc=G, burnin=G, thin=1, verbose=G,
    
51. b0=b0, B0=B0, c0=c0, d0=d0, r0=r0, R0=R0)
    
52. ## latent state changes
    
53. plotState(out1)
    
54. ## print mcmc output
    
55. summary(out1)

复制代码

https://cran.r-project.org/web/packages/MCMCpack/MCMCpack.pdf

好帖支持