混合模型学习笔记3

################################################################
### 02429 - Analysis of correlated data: Mixed Linear Models ###
### R-script for eNote-2                                     ###
################################################################

require(diagram) # For factor structure diagrams

## Create the list of factor names with indices
names <- c(expression("[I]"[24]^{30}),
           expression(atm:temp[2]^{6}),
           expression(atm[1]^{2}),
           expression(temp[2]^{3}),
           expression(0[1]^{1}))

## As there are 5 factors, create the 5x5 matrix of zeros
M <- matrix(nrow = 5, ncol = 5, byrow = TRUE, data = 0)

## Envision the structure: e.g., I need an arrow from the first
## factor in my list to the second, so assign something to M[2,1]
M[2, 1] <- M[3,2] <- M[4,2]  <- M[5,3] <- M[5,4] <- ""

## Make the diagram:
plotmat(M, pos = c(1, 1, 2, 1), name = names, lwd = 2,
        box.lwd = 1, cex.txt = 1, box.size = 0.1, box.type = "square",
        box.prop = 0.4, arr.type = "triangle", curve = 0)

## Function for rotation of factor structure diagram
matrix_position <- function(pos_vec) {
   n <- sum(pos_vec) # rows
   m <- length(pos_vec)-2 # inner layers
   d_hori <- 0.8/(m+1)
   bot <- 0.1; mid <- 0.5; top <- 0.9
   pos_mat <- matrix(nrow=n, ncol=2)
   pos_mat[1,1] <- bot; pos_mat[n,1] <- top
   pos_mat[1,2] <- pos_mat[n,2] <- mid
   cum_pos <- cumsum(pos_vec)
   for (i in 1:m) {
     n_vert <- pos_vec[i+1]
     d_vert <- 0.8 / (n_vert + 1)
     for (j in 1:n_vert){
       pos_mat[cum_pos[i] + j, 2] <- 0.1 + j*d_vert
       pos_mat[cum_pos[i] + j, 1] <- 0.1 + i*d_hori
     }
   }
   return(pos_mat)
}

## Rotated factor structure diagram
plotmat(M, pos = matrix_position(c(1, 1, 2, 1)), name = names, lwd = 2,
        box.lwd = 1, cex.txt = 1, box.size = 0.1, shadow.size = 0,
        box.type = "square", box.prop = 0.5, arr.type = "triangle",
        curve = 0)

## Create the list of factor names with indices
names <- c(expression("[I]"[266]^{300}),
           expression(depth:width[8]^{15}),
           expression("[plank]"[19]^{20}),
           expression(width[2]^{3}),
           expression(depth[4]^{5}),
           expression(0[1]^{1}))

## Since there are 6 factors, create the 6x6 matrix of zeros
M <- matrix(nrow = 6, ncol = 6, byrow = TRUE, data = 0)

## Envision the structure: E.g., I need an arrow from the first factor level to the second,
## so assign something to M[2,1], etc.
M[2, 1]  <- M[3, 1] <- M[4, 2]  <- M[5, 2] <- ""
M[6, 3] <- M[6, 4] <- M[6, 5] <- ""

## Make the diagram
plotmat(M, pos = c(1, 1, 3, 1), name = names, lwd = 2,
        box.lwd = 1, cex.txt = 1, box.size = 0.1, shadow.size = 0,
        box.type = "square", box.prop = 0.4, arr.type = "triangle",curve=0)

## Create the list of factor names with indices
names <- c(expression("[I]"[24]^{30}),
           expression(atm:temp[2]^{6}),
           expression(atm[1]^{2}),
           expression(temp[2]^{3}),
           expression(0[1]^{1}))

## Define the coordinates of all the terms in a grid, here (2x8; 2x8):
x <- c(2, 4, 6, 6, 8)
y <- c(5, 5, 3, 7, 5)

## Make an empty plot without any arrows or text:
plot(NA, NA, xlim = c(2, 8), ylim = c(2, 8), type = "n", axes = F,
     xlab = "", ylab = "")
text(x, y, names) # Add text according to the 'names' vector

## Define coordinates for the beginning and end points,
## (x0, y0) and (x1, y1), for 5 arrows:
x0 <- c(1.8, 4.2, 4.2, 6, 6) + .5
y0 <- c(5, 5, 5, 3, 7)
x1 <- c(2.7, 5, 5, 7.2, 7.2) + .5
y1 <- c(5, 3, 7, 5, 5)
arrows(x0, y0, x1, y1, length = 0.1) # 'length' applies to arrrow head

draw.arrows <- function(n) {
  ret <- locator(2*n)
  o <- (1:(2*n)) %% 2 == 1
  e <- (1:(2*n)) %% 2 == 0
  x0 <- ret$x[o]
  y0 <- ret$y[o]
  x1 <- ret$x[e]
  y1 <- ret$y[e]
  return(list(x0 = x0, y0 = y0, x1 = x1, y1 = y1))
}

plot(NA, NA, xlim = c(2, 8), ylim = c(2, 8), type = "n", axes = F,
    xlab = "", ylab = "")
text(x, y, names)
arr <- draw.arrows(5)

arrows(arr$x0, arr$y0, arr$x1, arr$y1, length = .1)

pdf("./mydiagram.pdf") # This is the path/name_of_figure
plot(NA, NA, xlim = c(2, 8), ylim = c(2, 8), type = "n", axes = F,
   xlab = "", ylab = "")
text(x, y, names)
arrows(arr$x0, arr$y0, arr$x1, arr$y1, length = .1)
# Close graphics device:
dev.off() # Only now, the figure is actually saved.