SEM2DSM/lavaan.prediction.function.R at master · angelini75/SEM2DSM · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
library(lavaan)


predict_lavaan <- function(fit, newdata = NULL){
  stopifnot(inherits(fit, "lavaan"))

  #Make sure we can use this
  if(!inspect(fit, "meanstructure")) stop("Need to supply meanstructure = TRUE in fit\n")
  if(is.null(newdata)){
    newdata <- data.frame(inspect(fit, "data"))
    names(newdata) <- lavNames(fit)
  }

  if(length(lavNames(fit, type="lv"))!=0) stop("Does not currently work with latent variables\n")

  #check for new data
  if(sum(!(lavNames(fit, type="ov.x") %in% names(newdata)))>0) stop("Not all exogenous variables supplied!")

  #Add some new columns to newdata
  newdata$Intercept <- 1
  newdata[lavNames(fit, "ov.nox")] <- 0


  mod_df <- data.frame(lhs = fit@ParTable$lhs,
                       op = fit@ParTable$op,
                       rhs = fit@ParTable$rhs,
                       exo = fit@ParTable$exo,
                       est = fit@ParTable$est,
                       se = fit@ParTable$se,
                       stringsAsFactors=FALSE)

  #Drop covariances
  mod_df <- mod_df[-which(mod_df$op=="~~"),]
  mod_df[which(mod_df$op=="~1"),]$rhs <- "Intercept"

  #get rid of exogenous on lhs
  mod_df <- mod_df[-which(mod_df$exo==1),]

  #Order by lhs
  mod_df <- mod_df[sort(mod_df$lhs, index.return=TRUE)$ix,]

  #let us know which variables on the rhs are exogenous
  mod_df$ord <- 0
  mod_df[which(!(mod_df$rhs %in% mod_df$lhs)),]$ord <- 1

  #Make a "order"
  ord_current <- 1
  while(sum(mod_df$ord==0)>0){
    for(r in unique(mod_df$lhs)){
      val <-  sum(mod_df[which(mod_df$lhs==r),]$ord==0)
      if(val==0) {
        mod_df[which(mod_df$lhs==r),]$ord <- ord_current
        mod_df[which(mod_df$rhs==r),]$ord <- ord_current+1
      }
    }
  ord_current <- ord_current +1
  }

  #correct for ragged ordering
  for(r in unique(mod_df$lhs)){
    mod_df[which(mod_df$lhs==r),]$ord <- max(mod_df[which(mod_df$lhs==r),]$ord)
  }

  #sort by order
  mod_df <- mod_df[sort(mod_df$ord, index.return=TRUE)$ix,]

  #now do the fitting in order
  fit_df <- data.frame(base = rep(1, nrow(newdata)))

  for(r in unique(mod_df$lhs)){
    subdf <- subset(mod_df, mod_df$lhs==r)
    #make a formula
    rhs <- paste0(subdf$rhs, collapse=" + ")
    form <- as.formula(paste0(r, " ~ ", rhs))

    #use formula to get right part of the data in right format
    mod_mat <- model.matrix(form, newdata)[,-1]
    new_val = mod_mat %*% subdf$est

    fit_df[[r]] <- new_val
    newdata[[r]] <- new_val
  }

  return(fit_df[,-1])

}


fitted_lavaan <- function(fit){
  predict_lavaan(fit)
}

residuals_lavaan <- function(fit){
  fitted_vals <- fitted_lavaan(fit)

  rawdata <- data.frame(inspect(fit, "data"))
  names(rawdata) <- lavNames(fit)

  res <- data.frame(base = rep(1, nrow(rawdata)))
  for(vals in names(fitted_vals)){
    res[[vals]] <- rawdata[[vals]] - fitted_vals[[vals]]
  }

  return(res[,-1])
}


########### EXAMPLE CODE

iris_mod <- "
Sepal.Width ~ Sepal.Length
Petal.Width ~ Sepal.Width + Petal.Length

"

iris_fit <- sem(iris_mod, data=iris, meanstructure = TRUE)

fitted_lavaan(iris_fit)

res <- residuals_lavaan(iris_fit)

par(mfrow=c(1,2))
for(i in 1:2) {
  qqnorm(res[,i], main=names(res)[i])
  qqline(res[,i])
}
par(mfrow=c(1,1))

newdata <- data.frame(Sepal.Length=1:10, Petal.Length=11:20)
pred <- predict_lavaan(iris_fit, newdata=newdata)

plot(Petal.Width ~ Sepal.Length, data=cbind(newdata, pred))

## MLR
sep.mod <- lm(Sepal.Width ~ Sepal.Length, iris)
pet.mod <- lm(Petal.Width ~ Sepal.Length + Petal.Length, iris)

res2 <- data.frame(Sepal.Width=residuals.lm(sep.mod),
                   Petal.Width=residuals.lm(pet.mod))
par(mfrow=c(1,2),pty="s")
for(i in 1:2) {
  qqnorm(res2[,i], main=names(res2)[i])
  qqline(res2[,i])
}

pred2 <- pred
pred2$Sepal.Width <- predict(sep.mod, newdata=newdata)
pred2$Petal.Width <- predict(pet.mod, newdata=newdata)

plot(Petal.Width ~ Sepal.Length, data=iris, col="blue", ylim = c(0,4), xlim=c(4,10))
points(Petal.Width ~ Sepal.Length, data=cbind(newdata, pred))
points(Petal.Width ~ Sepal.Length, data=cbind(newdata, pred2), col="red")