[Sem-additions-commits] r5 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 21 08:38:23 CEST 2010
Author: jebyrnes
Date: 2010-04-21 08:38:21 +0200 (Wed, 21 Apr 2010)
New Revision: 5
Added:
pkg/R/add.paths.R
pkg/R/anova.adjchisq.r
pkg/R/update.mod.r
pkg/R/update.model.r
pkg/man/add.paths.Rd
pkg/man/anova.adjchisq.Rd
pkg/man/update.mod.Rd
pkg/man/update.model.Rd
Modified:
pkg/ChangeLog
pkg/DESCRIPTION
pkg/GOALS
pkg/R/Ktrans.R
pkg/R/adf.wmat.R
pkg/R/browneChisq.R
pkg/R/bsboot.R
pkg/R/delta_matrix.R
pkg/R/helpers.R
pkg/R/ml.wmat.R
pkg/R/model.cov.R
pkg/R/robust_se.R
pkg/R/robust_summary.R
pkg/R/sbDiff.R
pkg/R/sbchisq.R
pkg/man/adjchisq.Rd
pkg/man/icmethod.Rd
pkg/man/plot.sem.Rd
pkg/man/rmsea.power.Rd
pkg/man/sbDiff.Rd
pkg/man/sbMeanChisq.Rd
pkg/man/sbchisq.Rd
pkg/man/sem.additions.Rd
Log:
Modified: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/ChangeLog 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,9 +1,23 @@
+2010-04-20 JEB
+ - added anova.adjchisq which uses the SB corrected ANOVA for nested model comparison
+ - now correct the forth moment by n instead of n-1 (suggestion from Yves Rosseel)
+ - fixed delta_matrix bug of not adjusting for correlated disturbance on both parts of P matrix
+ -added update.mod as a general way to update a ram object
+ -added methods to deal with lavaan models (sets of formulas)
+
+2010-02-09 JEB
+ - now using matrixcalc for creating of duplication matrix
+ - fixed bug in p_deriv_mat, no longer including fixed params
+ in calculation of delta matrix
+
2009-12-29 JEB
- fixed a mistaken generalized inverse in SB chisq
- fixed robust_summary to accomodate fixed parameters
+
2009-05-12 JEB
- Multiple fixes to Rd files.
- Added rmsea power test by Joerg Evermann
+
2008-09-02 JEB
- Error in delta_matrix fixed. Delta matrix now largely squares with STATA output from Stas Kolenikov
- Added the Chi Square Difference test
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/DESCRIPTION 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,10 +1,10 @@
Package: sem.additions
-Version: 0.1-05
-Date: 2009-1229
+Version: 0.1-10
+Date: 2010-04-20
Title: Additional Methods for Structural Equation Modelling
Author: Jarrett Byrnes <byrnes at msi.ucsb.edu>
Maintainer: Jarrett Byrnes <byrnes at msi.ucsb.edu>
-Depends: R(>= 2.7), sem, mvnormtest, sna, boot
+Depends: R(>= 2.8), sem, mvnormtest, sna, boot, matrixcalc, MASS, igraph
Description: A series of methods for easier manipulation of RAM model objects, graphing tools for sem objects, and a variety of fit indices and corrections for non-normality.
License: GPL-3
URL: https://r-forge.r-project.org/projects/sem-additions/
Modified: pkg/GOALS
===================================================================
--- pkg/GOALS 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/GOALS 2010-04-21 06:38:21 UTC (rev 5)
@@ -3,10 +3,12 @@
1.0 - Code would include:
- Model manipulation routines
- SB Chisq, SB Mean Chi Square, YB Chi Square, BCC, ECVI, and more
- SB Robust Standard Errors and Sandwich Standard Errors
- SB ANOVA
- Multigroup Analysis
+ Model manipulation routines - done
+ SB Chisq, SB Mean Chi Square - done
+ YB Chi Square, BCC, ECVI, and more - many done
+ SB Robust Standard Errors - done
+ Sandwich Standard Errors
+ SB ANOVA - done
+ Multigroup Analysis - done for equal sample size, but not yet included
Graphical Interface for Model Specification
Retrospective power analysis
Modified: pkg/R/Ktrans.R
===================================================================
--- pkg/R/Ktrans.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/Ktrans.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -4,15 +4,18 @@
# symmetric matrix & placing them in a vector
# Adapted from Stas Kolenikov's MATA code for a duplication matrix by JEB
#
-# Last Updated 9/02/08
+# 2/8/20 - realized matrixcalc package has a MUCH faster way of doing this
+# Last Updated 2/08/10
##
Ktrans<-function(num.vars){
- pstar<-num.vars*(num.vars+1)/2
- Ipstar<-ident.mat(pstar)
- D.mat<-c()
- for (k in 1:pstar) D.mat<-c(D.mat, vec(invvech(Ipstar[,k])))
- D.mat<-matrix(D.mat, ncol=pstar)
+ #pstar<-num.vars*(num.vars+1)/2
+ #Ipstar<-ident.mat(pstar)
+ #D.mat<-c()
+ #for (k in 1:pstar) D.mat<-c(D.mat, vec(invvech(Ipstar[,k])))
+ #D.mat<-matrix(D.mat, ncol=pstar)
+ require(matrixcalc)
+ D.mat<-duplication.matrix(num.vars)
return(D.mat)
}
Added: pkg/R/add.paths.R
===================================================================
--- pkg/R/add.paths.R (rev 0)
+++ pkg/R/add.paths.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -0,0 +1,15 @@
+add.paths<-function(from, to, values=rep(NA, length(from)*length(to))){
+ newram<-matrix(rep(NA, 3*length(from)*length(to)), ncol=3)
+ par.index<-0
+ for (i in 1:length(to)){
+ for (j in 1:length(from)){
+ par.index<-par.index+1
+ newram[par.index,1]<-paste(from[j], "->", to[i], sep="")
+ newram[par.index,2]<-paste(from[j], ".", to[i], sep="")
+ newram[par.index,3]<-values[par.index]
+ }
+ }
+
+ class(newram)<-"mod"
+ return(newram)
+}
Modified: pkg/R/adf.wmat.R
===================================================================
--- pkg/R/adf.wmat.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/adf.wmat.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,6 +1,8 @@
#calculate ADF weight matrix from raw data
#this is actually waaaay faster, and produces correct output. The Bentler adapted code was giving me the wrong matrix for some reason...
+#last updated 4/20/10
+
adf.wmat<-function(raw_data){
n<-length(raw_data[,1])
n.col<-length(names(raw_data))
@@ -22,7 +24,7 @@
sc[[a.name]]<-z[,i]*z[,j]
}
- adf_mat<-var(data.frame(sc))
+ adf_mat<-var(data.frame(sc))* (n-1)/n #to correct for n-1
return(adf_mat)
}
Added: pkg/R/anova.adjchisq.r
===================================================================
--- pkg/R/anova.adjchisq.r (rev 0)
+++ pkg/R/anova.adjchisq.r 2010-04-21 06:38:21 UTC (rev 5)
@@ -0,0 +1,49 @@
+###
+#from http://www.statmodel.com/chidiff.shtml
+# Satorra-bentler adjusted chi sq
+# last updated 4/20/10
+###
+anova.adjchisq<-function(adjobj0, adjobj1){
+
+ #create an output table
+ out.mat<-matrix(ncol=5, nrow=2)
+ rownames(out.mat)<-c("Model 1", "Model 2")
+ colnames(out.mat)<-c("Model DF", "Model Chisq", "Df", "LR Chisq", "Pr(>Chisq)")
+
+ out.mat[1,1]<-adjobj0$df
+ out.mat[1,2]<-adjobj0$chisq.scaled
+
+ out.mat[2,1]<-adjobj1$df
+ out.mat[2,2]<-adjobj1$chisq.scaled
+
+ #switching to get order right
+ sbs.nested<-adjobj0
+ sbs.full<-adjobj1
+ if(adjobj0$df<adjobj1$df){
+ sbs.nested<-adjobj1
+ sbs.full<-adjobj0
+ }
+
+
+
+ t0<-sbs.nested$chisq
+ tr0<-sbs.nested$chisq.scaled
+ t1<-sbs.full$chisq
+ tr1<-sbs.full$chisq.scaled
+
+ c0<-sbs.nested$c
+ c1<-sbs.full$c
+
+ d0<-sbs.nested$df
+ d1<-sbs.full$df
+
+ cd = (d0 * c0 - d1*c1)/(d0 - d1)
+ trd = (t0 - t1)/cd
+
+ out.mat[2,3]<-d0-d1
+ out.mat[2,4]<-trd
+ out.mat[2,5]<-1-pchisq(trd, d0-d1)
+
+ print(out.mat, na.print = "")
+
+ }
\ No newline at end of file
Modified: pkg/R/browneChisq.R
===================================================================
--- pkg/R/browneChisq.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/browneChisq.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,3 +1,5 @@
+#last updated 4/20/2010
+
browneChisq<-function(sem.obj, sem.data, ...){
adj.obj<-sbchisq(sem.obj, sem.data, ...)
adj.obj$w_mat<-NA
@@ -12,4 +14,4 @@
adj.obj$chisq.scaled<- adj.obj$N * t(s-sigma) %*% adj.obj$res_u %*% (s-sigma)
adj.obj$p<-1-pchisq(adj.obj$chisq.scaled, adj.obj$df)
return(adj.obj)
-}
+}
\ No newline at end of file
Modified: pkg/R/bsboot.R
===================================================================
--- pkg/R/bsboot.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/bsboot.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -98,7 +98,7 @@
summary.bsboot<-function(object){
mean.chisq<-mean(object$chisq.vec, na.rm=T)
- se.chisq<-se(object$chisq.vec)
+ se.chisq<-sd(object$chisq.vec, na.rm=T)/sqrt(length(na.omit(object$chisq.vec)))
var.chisq<-var(object$chisq.vec, na.rm=T)
p<-1-pchisq(mean.chisq, object$df)
ret.mat<-matrix(c(mean.chisq, var.chisq, se.chisq, object$R, object$df, p), nrow=1)
Modified: pkg/R/delta_matrix.R
===================================================================
--- pkg/R/delta_matrix.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/delta_matrix.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -3,13 +3,18 @@
# respect to the non-duplicated fitted covariances
# Adapted from SPSS code by Peter Bentler
#
-# Changelog
+# Changelog
+# 3/5/10 - fixed bug of not adjusting for correlated disturbance on both parts of P matrix
+# 2/8/10 - fixed bug of using fixed params
# 9/1/08 - fixed matrix selection bug
##
delta_matrix<-function(sem.object, adj=1e-04){
p.star<-sem.object$n*(sem.object$n+1)/2
- nparams<-length(sem.object$ram[,1])
+ #nparams<-length(sem.object$ram[,1])
+ #only use free parameters
+ par.idx<-which(sem.object$ram[,"parameter"]>0)
+ nparams<-length(par.idx)
delta.mat<-matrix(0,nparams, p.star)
rownames(delta.mat)<-rep(NA, nparams)
colnames(delta.mat)<-vech(matrix.names(sem.object$C))
@@ -19,20 +24,27 @@
J<-sem.object$J
m<-sem.object$m
- #iterate through the ram
- for (i in 1:nparams){
- A<-sem.object$A
- P<-sem.object$P
- if(sem.object$ram[i,1]==1){
- A[sem.object$ram[i,2],sem.object$ram[i,3]]<-
- A[sem.object$ram[i,2],sem.object$ram[i,3]]+adj
- }else{
- P[sem.object$ram[i,2],sem.object$ram[i,3]]<-
- P[sem.object$ram[i,2],sem.object$ram[i,3]]+adj
- }
+ #iterate through the ram
+ #but only look at free params
+ for (i in 1:nparams){
+ A<-sem.object$A
+ P<-sem.object$P
+ from<-sem.object$ram[par.idx[i],2]
+ to<-sem.object$ram[par.idx[i],3]
+ path_type <-sem.object$ram[par.idx[i],1]
+ if(path_type==1){
+ A[from, to] <-
+ A[from, to] + adj
+ }else{
+ P[from, to]<- P[from,to]+adj
+
+
+ #symmetric P matrix
+ P[to, from]<- P[to,from]
+ }
+
+ I.Ainv <- solve(diag(m) - A)
- I.Ainv <- solve(diag(m) - A)
-
#calculate fitted covarianve matrix from RAM form
#C= J(Im A) 1P[(Im A) 1]J using the RAM formulation
C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J)
@@ -41,11 +53,11 @@
#get difference with model cov matrix / adj
#linearize, and put into delta matrix
delta.mat[i,]<-(C.vech - C.vect)/adj
- rownames(delta.mat)[i] <-rownames(sem.object$ram)[i]
+ rownames(delta.mat)[i] <-rownames(sem.object$ram)[par.idx[i]]
if(rownames(delta.mat)[i] ==""){
rownames(delta.mat)[i] <-
- paste(vars[sem.object$ram[i,3]],
- vars[sem.object$ram[i,2]],
+ paste(vars[sem.object$ram[par.idx[i],3]],
+ vars[sem.object$ram[par.idx[i],2]],
sep="-")
}
Modified: pkg/R/helpers.R
===================================================================
--- pkg/R/helpers.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/helpers.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -26,13 +26,14 @@
#vech, as the ks vech borks on certain valid matrices due to type issues
# code based on ks library
-vech<-function (x)
-{
- d <- ncol(x)
- vechx <- vector()
- for (j in 1:d) vechx <- c(vechx, x[j:d, j])
- return(vechx)
-}
+#now using matrixcalc code
+#vech<-function (x)
+#{
+# d <- ncol(x)
+# vechx <- vector()
+# for (j in 1:d) vechx <- c(vechx, x[j:d, j])
+# return(vechx)
+#}
#invvech, so as to remove dependency on ks - code based on ks library
invvech<-function (x)
@@ -45,15 +46,16 @@
return(invvechx)
}
-vec<-function (x, byrow = FALSE)
-{
- if (byrow)
- x <- t(x)
- d <- ncol(x)
- vecx <- vector()
- for (j in 1:d) vecx <- c(vecx, x[, j])
- return(vecx)
-}
+#now using code from matrixcalc
+#vec<-function (x, byrow = FALSE)
+#{
+# if (byrow)
+# x <- t(x)
+# d <- ncol(x)
+# vecx <- vector()
+# for (j in 1:d) vecx <- c(vecx, x[, j])
+# return(vecx)
+#}
#gets names of each row_col of a matrix
Modified: pkg/R/ml.wmat.R
===================================================================
--- pkg/R/ml.wmat.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/ml.wmat.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -14,7 +14,7 @@
if(useFit){An<-sem.obj$C} else {An<-sem.obj$S}
Dp<-Ktrans(p)
An.inv<-solve(An)
- w_mat<-2^-1 * t(Dp) %*% kronecker(An.inv, An.inv) %*% Dp
+ w_mat<-0.5 * t(Dp) %*% kronecker(An.inv, An.inv) %*% Dp
rownames(w_mat)<-vech(matrix.names(sem.obj$C))
colnames(w_mat)<-rownames(w_mat)
return(w_mat)
Modified: pkg/R/model.cov.R
===================================================================
--- pkg/R/model.cov.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/model.cov.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,7 +1,10 @@
-model.cov<-function(ram.object, input.data){
+#last modified 4/10/10
+#changelog
+# 4/10/10 - added na.rm argument
+model.cov<-function(ram.object, input.data, na.rm=T, use="complete.obs"){
newdata<-model.data(ram.object, input.data)
- cov.matrix<-var(newdata, na.rm=T)
+ cov.matrix<-var(newdata, na.rm=na.rm, use=use)
return(cov.matrix)
}
Modified: pkg/R/robust_se.R
===================================================================
--- pkg/R/robust_se.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/robust_se.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -8,21 +8,31 @@
return(ret)
}
-robust_se<-function(sem.obj, adj.obj=NA, data.obj=NA, useFit=F){
- if(is.na(adj.obj[1]) && is.na(data.obj)) stop ("Need a data or sbchisq object")
-
- #get all of those matrices we'll need
- if(is.na(adj.obj[1])){ adj.obj<-sbchisq(sem.obj, data.obj, useFit=useFit)}
-
- ncases<-sem.obj$N
-
- #calculate the hessian
- hes<-sem_hessian(adj.obj$w_mat, adj.obj$p_deriv_mat)
- info_m<-solve(hes)
-
- acov<- info_m %*% ( t(adj.obj$p_deriv_mat) %*% adj.obj$w_mat %*% adj.obj$w_adf %*% adj.obj$w_mat %*% adj.obj$p_deriv_mat) %*% info_m
-
- ret<-sqrt(diag(acov) / (ncases-1))
- names(ret)<-colnames(adj.obj$p_deriv_mat)
- return(ret)
+robust_se<-function (sem.obj, adj.obj = NA, data.obj = NA, useFit = F, useGinv=F)
+{
+ if (is.na(adj.obj[1]) && is.na(data.obj))
+ stop("Need a data or sbchisq object")
+ if (is.na(adj.obj[1])) {
+ adj.obj <- sbchisq(sem.obj, data.obj, useFit = useFit)
+ }
+ ncases <- sem.obj$N
+ hes <- sem_hessian(adj.obj$w_mat, adj.obj$p_deriv_mat)
+ #having issues with this inversion in certain models
+ #until it is locked down, ginv provides a reasonable
+ #approximation
+ info_m <- try(solve(hes), silent=T)
+ if(class(info_m)=="try-error" & useGinv==T) {
+ info_m <-ginv(hes)
+ ginvFlag<-T
+ }
+ acov <- info_m %*% (t(adj.obj$p_deriv_mat) %*% adj.obj$w_mat %*%
+ adj.obj$w_adf %*% adj.obj$w_mat %*% adj.obj$p_deriv_mat) %*%
+ info_m
+ ret <- sqrt(diag(acov)/(ncases - 1))
+ names(ret) <- colnames(adj.obj$p_deriv_mat) #not sure if this is needed
+ return(ret)
+
+ #have dealt with this in delta_mat
+ #return(ret[-which(sem.obj$ram[,4]==0)])
}
+
Modified: pkg/R/robust_summary.R
===================================================================
--- pkg/R/robust_summary.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/robust_summary.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,11 +1,12 @@
-robust_summary<-function (sem.obj, adj.obj = NA, data.obj = NA, useFit = F)
+#last updated 4/20/2010
+robust_summary<-function (sem.obj, adj.obj = NA, data.obj = NA, useFit = F, useGinv=F)
{
if (is.na(adj.obj[1]) && is.na(data.obj))
stop("Need a data or adjchisq object")
if (is.na(adj.obj[1])) {
- adj.obj <- sbchisq(sem.obj, data.obj, useFit = useFit)
+ adj.obj <- sbchisq(sem.obj, data.obj, useFit = useFit, useGinv= useGinv)
}
- ses <- robust_se(sem.obj, adj.obj = adj.obj)
+ ses <- robust_se(sem.obj, adj.obj = adj.obj, useGinv=adj.obj$ginvFlag)
se <- rep(NA, length(ses))
index<-0
for (i in 1:length(sem.obj$ram[, 1])) {
@@ -27,18 +28,3 @@
print(coef.mat)
}
-add.paths<-function(from, to, values=rep(NA, length(from)*length(to))){
- newram<-matrix(rep(NA, 3*length(from)*length(to)), ncol=3)
- par.index<-0
- for (i in 1:length(to)){
- for (j in 1:length(from)){
- par.index<-par.index+1
- newram[par.index,1]<-paste(from[j], "->", to[i], sep="")
- newram[par.index,2]<-paste(from[j], ".", to[i], sep="")
- newram[par.index,3]<-values[par.index]
- }
- }
-
- class(newram)<-"mod"
- return(newram)
-}
\ No newline at end of file
Modified: pkg/R/sbDiff.R
===================================================================
--- pkg/R/sbDiff.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/sbDiff.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -1,23 +1,25 @@
-#last modified 9/2/2008
+#last modified 4/20/2010
-sbDiff<-function(sbobj.A, sbobj.B){
-
- #guard against specifying the wrong order of who is nested
- if(sbobj.A$df < sbobj.B$df){m0<-sbobj.B; m1<-sbobj.A #B is nested
- } else {m0<-sbobj.B; m1<-sbobj.A} #A is nested
+sbDiff<-function(sbobj.A, sbobj.B) anova.adjchisq(sbobj.A, sbobj.B)
- df<-m0$df - m1$df
-
- #calculation of cd from http://www.statmodel.com/chidiff.shtml
- cd<-(m0$df * m0$c - m1$df * m1$c)/(m0$df - m1$df)
-
- TRd <- (m0$chisq - m1$chisq)/cd
-
- p<-1-pchisq(TRd, df)
-
- retmat<-matrix(c(TRd, df, p), nrow=1)
- rownames(retmat)<-""
- colnames(retmat)<-c("Chi Sq", "DF", "p")
-
- retmat
- }
\ No newline at end of file
+#sbDiff<-function(sbobj.A, sbobj.B){
+#
+# #guard against specifying the wrong order of who is nested
+# if(sbobj.A$df < sbobj.B$df){m0<-sbobj.B; m1<-sbobj.A #B is nested
+# } else {m0<-sbobj.B; m1<-sbobj.A} #A is nested
+#
+# df<-m0$df - m1$df
+#
+# #calculation of cd from http://www.statmodel.com/chidiff.shtml
+# cd<-(m0$df * m0$c - m1$df * m1$c)/(m0$df - m1$df)
+#
+# TRd <- (m0$chisq - m1$chisq)/cd
+#
+# p<-1-pchisq(TRd, df)
+#
+# retmat<-matrix(c(TRd, df, p), nrow=1)
+# rownames(retmat)<-""
+# colnames(retmat)<-c("Chi Sq", "DF", "p")
+#
+# retmat
+# }
\ No newline at end of file
Modified: pkg/R/sbchisq.R
===================================================================
--- pkg/R/sbchisq.R 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/R/sbchisq.R 2010-04-21 06:38:21 UTC (rev 5)
@@ -2,47 +2,39 @@
#Calcute the Satorra-Bentler Chi Square Index
#Adapted from SPSS code by Peter Bentler
#Returns a number of matrices also useful to other functions
+#last updated 4/20/2010
##
-sbchisq<-function(sem.obj, sem.data, adj=1e-04, useFit=F){
- props<-sem.props(sem.obj)
- sem.prop <- props$chisq
- chisq<-props$chisq
- df<-props$df
-
- #calculate the ADF weight matrix
- w_adf<-adf.wmat(sem.data)
-
- #get the GLS weight matrix
- w_mat<-ml.wmat(sem.obj, useFit=useFit)
-
- #get the delta matrix, the partial derivative of changes in covariances
- #for changes in parameter values
- p_deriv_mat<- delta_matrix(sem.obj)
-
- #calculate the LS residual weight matrix
- res_u<- w_mat - (w_mat%*%p_deriv_mat %*% solve(t(p_deriv_mat) %*% w_mat %*% p_deriv_mat) %*% t(p_deriv_mat) %*% w_mat)
-
- #compute scaling statistic
- ug<-res_u %*% w_adf
- scale_stat<-tr(ug)/df
-
- #compute
- chisq.scaled<-chisq/scale_stat
- p.old<-1-pchisq(chisq, df)
- p<-1-pchisq(chisq.scaled, df)
- ret<-list( chisq = chisq,
- t = sem.obj$t,
- df = df,
- p.old = p.old,
- c = scale_stat,
- chisq.scaled = chisq.scaled,
- p = p,
- w_mat = w_mat,
- p_deriv_mat = p_deriv_mat,
- w_adf = w_adf,
- res_u = res_u,
- N = length(sem.data[,1]))
- class(ret)<-"adjchisq"
-
- return(ret)
+sbchisq<-function (sem.obj, sem.data, adj = 1e-04, useFit = F, useGinv=F)
+{
+ props <- sem.props(sem.obj)
+ sem.prop <- props$chisq
+ chisq <- props$chisq
+ df <- props$df
+ w_adf <- adf.wmat(sem.data)
+ w_mat <- ml.wmat(sem.obj, useFit = useFit)
+ p_deriv_mat <- delta_matrix(sem.obj)
+
+ #having issues with this inversion in certain models
+ #until it is locked down, ginv provides a reasonable
+ #approximation
+ ginvFlag<-F
+ invMat<-try(solve(t(p_deriv_mat) %*% w_mat %*% p_deriv_mat), silent=T)
+ if(class(invMat)=="try-error" & useGinv==T) {
+ invMat<-ginv(t(p_deriv_mat) %*% w_mat %*% p_deriv_mat)
+ ginvFlag<-T
+ }
+ res_u <- w_mat - (w_mat %*% p_deriv_mat %*% invMat %*% t(p_deriv_mat) %*% w_mat)
+
+ ug <- res_u %*% w_adf
+ scale_stat<-sum(diag(ug))/df
+ chisq.scaled <- chisq/scale_stat
+ p.old <- 1 - pchisq(chisq, df)
+ p <- 1 - pchisq(chisq.scaled, df)
+
+ ret <- list(chisq = chisq, t = sem.obj$t, df = df, p.old = p.old,
+ c = scale_stat, chisq.scaled = chisq.scaled, p = p, w_mat = w_mat,
+ p_deriv_mat = p_deriv_mat, w_adf = w_adf, res_u = res_u,
+ N = length(sem.data[, 1]), ginvFlag=ginvFlag)
+ class(ret) <- "adjchisq"
+ return(ret)
}
Added: pkg/R/update.mod.r
===================================================================
--- pkg/R/update.mod.r (rev 0)
+++ pkg/R/update.mod.r 2010-04-21 06:38:21 UTC (rev 5)
@@ -0,0 +1,122 @@
+#####
+# update.mod
+# updates a ram model
+# by addin, deleting, or replacing elements
+#
+# uses scan to get a list of modifications
+#
+# to delete, the format for a line is
+# delete, NAME, SORT
+# where name is the name of a variable, path, or coefficient
+# and sort is one of those three words
+#
+# to add a new path, simply use the following format
+# add, a->b, coef.name, start.value
+# where, aside from the word add, the rest is the same as a line
+# in specify.model()
+#
+# to replace one variable with a different one, use the following
+# replace, OLD, NEW
+# where old is the name of a variable currently in the model
+# and new is the name of the variable to which you would like to change it
+#
+# Last updated 3/1/2010
+#####
+
+update.mod<-function(object, file=""){
+ modmat<-scan(file = file, what = list(change="", var1 = "", var2 = "",
+ var3 = "", var4=""), sep = ",", strip.white = TRUE,
+ comment.char = "#", fill = TRUE)
+ modmat<-cbind(modmat$change, modmat$var1, modmat$var2, modmat$var3)
+
+ #####ADDITIONS
+
+ if ("add" %in% modmat[,1]){
+ addmat<-modmat[which(modmat[,1]=="add"), 2:4]
+
+ #length check in case only one thing is being added
+ if(length(addmat)==3) addmat <-matrix(addmat, ncol=3)
+ class(addmat)<-"mod"
+
+ #combine the models
+ object<-combine.models(object, addmat)
+ }
+
+ #######DELETIONS
+
+ ######
+ # delete.model.element
+ #
+ # a local function
+ # takes care of deleting elements from a ram object
+ ######
+
+ delete.model.element<- function (delete.text, old.model, type = "path") {
+ #a failed match provides an infomative enough error message
+ type<-match.arg(type, c("path", "variable", "coefficient"))
+ col.index<-list("path"=1, "variable"=1, "coefficient"=2)[[type]]
+ delete.text <- strip.white(delete.text)
+ old.model<-old.model[which(is.na(pmatch(strip.white(old.model[,col.index]), delete.text))),]
+ class(old.model)<-"mod"
+ return(old.model)
+ }
+
+ if("delete" %in% modmat[,1]){
+
+ deletemat<-modmat[which(modmat[,1]=="delete"), 2:3]
+ #length check in case only one thing is being deleted
+ if(length(deletemat)==2) deletemat <-matrix(deletemat, ncol=2)
+
+ #scroll through and delete all relevant entities
+ for (i in 1:length(deletemat[,1]))
+ object<-delete.model.element(deletemat[i,1], object, deletemat[i,2])
+ }
+ #delete.model.element("theta1", model.dhp, type="variable")
+
+ #############REPLACE
+ #to substitude one variable in for another
+
+ if("replace" %in% modmat[,1]){
+ submat<-modmat[which(modmat[,1]=="replace"),2:3]
+ #the variable is just to keep things quiet
+ mapply( function(x,y) object[,1:2]<<-gsub(x, y, object[,1:2]), submat[,1], submat[,2])
+ }
+ return(object)
+}
+
+
+
+#a quick demo
+#model.dhp <- specify.model()
+# RParAsp -> RGenAsp, gam11, NA
+# RIQ -> RGenAsp, gam12, NA
+# RSES -> RGenAsp, gam13, NA
+# FSES -> RGenAsp, gam14, NA
+# RSES -> FGenAsp, gam23, NA
+# FSES -> FGenAsp, gam24, NA
+# FIQ -> FGenAsp, gam25, NA
+# FParAsp -> FGenAsp, gam26, NA
+# FGenAsp -> RGenAsp, beta12, NA
+# RGenAsp -> FGenAsp, beta21, NA
+# RGenAsp -> ROccAsp, NA, 1
+# RGenAsp -> REdAsp, lam21, NA
+# FGenAsp -> FOccAsp, NA, 1
+# FGenAsp -> FEdAsp, lam42, NA
+# RGenAsp <-> RGenAsp, ps11, NA
+# FGenAsp <-> FGenAsp, ps22, NA
+# RGenAsp <-> FGenAsp, ps12, NA
+# ROccAsp <-> ROccAsp, theta1, NA
+# REdAsp <-> REdAsp, theta2, NA
+# FOccAsp <-> FOccAsp, theta3, NA
+# FEdAsp <-> FEdAsp, theta4, NA
+#
+#update(model.dhp)
+#add, newvar -> FEdASP, gam71, NA
+#add, newvar <-> newvar, theta5, NA
+#delete, FSES, variable
+#delete, theta1, variable #this should fail
+#delete, theta2, coefficient #this should not
+#delete, RIQ -> RGenAsp, path
+#replace, FEdAsp, HELLO!
+#replace, ROccAsp, GOODBYE
+
Added: pkg/R/update.model.r
===================================================================
--- pkg/R/update.model.r (rev 0)
+++ pkg/R/update.model.r 2010-04-21 06:38:21 UTC (rev 5)
@@ -0,0 +1,213 @@
+
+#takes a model and splits it into lines
+split.model<-function(a.model){
+ modvec<-strsplit(a.model, "\n")[[1]]
+ return(modvec)
+}
+
+#clean up a model line and return a formula only
+clean.modline<-function(modline){
+ modline<-gsub("=", "", modline)
+ modline<-gsub("~~", "~", modline)
+ return(as.formula(modline))
+}
+
+#get the response variable of a formula
+response<-function(a.formula){
+ gsub("()","", a.formula[2])
+ }
+
+#get the predictor variables of a formula
+predictors<-function(a.formula){
+# terms(a.formula)[[3]]
+ gsub("()","", a.formula[3])
+ }
+
+#what kind of operator is a formula using?
+getOperator<-function(a.formula){
+ op<-NA
+ for (a.op in c("~", "~~", "=~"))
+ if(grepl(a.op, a.formula, perl=T)) op<-a.op
+ return(op)
+ }
+
+
+#####
+#update a lavaan model
+#####
+update.model<-function(old.model, changes){
+ #first, break changes into list of response variables, operators, and updates
+ #then put it into a list for comparison against the old model
+ changes<-split.model(changes)
+ changelist<-list()
+ for (i in 1:length(changes)){
+ op<-getOperator(changes[i])
+ if(!is.na(op)){
+ a.formula<-clean.modline(changes[i])
+ changelist[[paste(response(a.formula), op, sep="")]]<-predictors(a.formula)
+ }
+ }
+
+ #now, run through the old model, and change things as needed
+ old.model<-split.model(old.model)
+ newmodel<-""
+ for(i in 1:length(old.model)){
+ #first, see if anything in the list matches
+ newline<-old.model[i]
+ op<-getOperator(old.model[i])
+ if(!is.na(op)){
+ a.formula<-clean.modline(old.model[i])
+ lookup<-changelist[[paste(response(a.formula), op, sep="")]]
+ if(!is.null(lookup)){
+ #OK, now, update the formula
+ a.formula<-update(a.formula, paste(response(a.formula),"~", lookup, sep=""))
+
+ #now make a new line
+ #strips off comments and other whitespacing...fix later
+ newline<-paste(response(a.formula), op, predictors(a.formula))
+
+ #making sure we have't stripped out everything
+ #this will mess with means only
+ #fix this later
+ if(predictors(a.formula)==1) newline=""
+
+ }
+
+ }
+ newmodel<-paste(newmodel, newline, "\n", sep="")
+ }
+
+ return(newmodel)
+}
+#
+###example
+#a.model<-'y~x+z
+# x~z
+# z=~a+b+c
+# y~~c
+#'
+#
+#some.changes<-'y~.-x
+# y~~.-c
+# x~.+y
+#'
+#
+#update.model(a.model, some.changes)
+
+#####
+#turn a lavaan model
+#into a ram object for the sem package
+#####
+model.to.ram<-function(a.model){
+ require(sem.additions)
+ a.model<-split.model(a.model)
+ ram<-matrix(ncol=3)
+
+ #for each line in the model
+ for(i in 1:length(a.model)){
+ op<-getOperator(a.model[i])
+
+ #make sure this is really a model line
+ #and now a comment or whitespace
+ if(!is.na(op)){
+
+ #what kind of path is this?
+ pathType<-"->"
+ if(op=="~~") pathType <-"<->"
+ a.formula<-clean.modline(a.model[i])
+ y<-response(a.formula)
+ x<-strip.white(strsplit(predictors(a.formula), "+", fixed=T)[[1]])
+
+ #iterate over each predictor variable
+ for(j in 1:length(x)){
+ #construct a new path
+ path<-paste(x[j], y, sep=pathType)
+
+ #see if this is a fixed parameter values
+ isFixed<-strsplit(x[j], '*', fixed=T)[[1]]
+ if(length(isFixed)>1) {
+ coef<-NA
+ startValue<-isFixed[1]
+ #otherwise....
+ }else{
+ #create a new coefficient name
+ coef<-paste(x[j], y, sep=".")
+ if(pathType=="<->") coef<-paste(coef, "cov", sep=".")
+ startValue<-NA
+ }
+ ram<-rbind(ram, c(path, coef, startValue))
+ }
+ }
+ }
+ #clean off the NA row
+ ram<-ram[2:length(ram[,1]),]
+ class(ram)<-"mod"
+ #add disturbance terms where needed
+ #using sem.additions code
+ ram<-add.errors(ram)
+ return(ram)
+}
+
+###example
+#a.model<-'y~x+z
+# x~3*z
+# z=~a+b+c
+# y~~c
+#'
+#
+#cat(a.model)
+#model.to.ram(a.model)
+
+
+###
+#can we go back?
+#we'll need to know which are latent variables
+#as there is no way of telling from the ram formulation
+###
+ram.to.model<-function(ram.obj, latents=NA){
+ # a function to determine if a value is fixed
+ #and return the proper predictor for a model
+ getFix<-function(a.row){
+ if(is.na(a.row[2])) return(paste(a.row[3], "*", sep=""))
+ return("")
+ }
+
+ model_list<-list()
+ for(i in 1:length(ram.obj[,1])){
+ ram.obj[i,1]<-strip.white(ram.obj[i,1])
+ iscov<-grepl("<->", ram.obj[i,1], perl=T)
+ vars<-strsplit(ram.obj[i,1], c("->", "<->")[as.numeric(iscov)+1])[[1]]
+
+ op<-"~"
+ #is this a latent variable
+ if(vars[2] %in% latents) op<-"=~"
+
+ #is this a covariance term
+ if(iscov) op<-"~~"
+
+ newKey<-paste(vars[2], op, sep="")
+
+ #is it in the model list?
+ if(is.null(model_list[[newKey]])){
+ #if not, start a new entry in the model lists
+ model_list[[newKey]]<-paste(getFix(ram.obj[i,]), vars[1], sep="")
+ }else{
+ #add the variable otherwise
+ model_list[[newKey]]<-paste(model_list[[newKey]], "+",getFix(ram.obj[i,]), vars[1], sep="")
+ }
+ }
+ a.mod<-""
+ for (j in 1:length(names(model_list))){
+ a.mod<-paste(a.mod, names(model_list)[j], model_list[[j]], "\n", sep="")
+ }
+ return(a.mod)
+}
+
+###example
+#a.ram<-specify.model()
+#a->b, ab, NA
+#b->c, bc, NA
+#d->c, NA, 3
+#a<->c, ac.cov, NA
+#
+#cat(ram.to.model(a.ram, latents="b"))
\ No newline at end of file
Added: pkg/man/add.paths.Rd
===================================================================
--- pkg/man/add.paths.Rd (rev 0)
+++ pkg/man/add.paths.Rd 2010-04-21 06:38:21 UTC (rev 5)
@@ -0,0 +1,45 @@
+\name{add.paths}
+\alias{add.paths}
+\title{Adds Paths From One Group of Variables to Another}
+\description{ Automates the adding of directional paths to ram models, rather than coding it all by hand. }
+
+\usage{add.paths(from, to, values=rep(NA, length(from)*length(to)))}
+\arguments{
+ \item{from}{A vector of predictor variables.}
+ \item{to}{A vector of endogenous response variables.}
+ \item{values}{Start values for the paths}
+}
+
+\details{Creates paths for all combinations of variables in one vector connecting to another. Useful for large models with many connections from certain predictor variables.}
+
+
+\seealso{
+ \code{\link{sem}}
+ \code{\link{specify.model}}
+}
+
+\examples{
+
+#a piece of a larger model
+
+
+model.kerch <- specify.model()
+ Intelligence -> Grades, gam51, NA
+ Siblings -> Grades, gam52, NA
+ FatherEd -> Grades, gam53, NA
+ FatherOcc -> Grades, gam54, NA
+ Intelligence -> EducExp, gam61, NA
+ Siblings -> EducExp, gam62, NA
+ FatherEd -> EducExp, gam63, NA
+ FatherOcc -> EducExp, gam64, NA
+
+
+model.kerch2<-add.paths(c("Intelligence", "Siblings", "FatherEd", "FatherOcc"),
+ c("Grades", "EducExp"))
+
+##compare the two submodels - they are the same
+##although with different names for coefficients
+model.kerch
+model.kerch2
+
+}
\ No newline at end of file
Modified: pkg/man/adjchisq.Rd
===================================================================
--- pkg/man/adjchisq.Rd 2010-02-02 19:59:42 UTC (rev 4)
+++ pkg/man/adjchisq.Rd 2010-04-21 06:38:21 UTC (rev 5)
@@ -5,7 +5,11 @@
\title{Adjusted Chi Squared Tests for Structural Equation Models}
\description{ Adjusted Chi Squared Tests for Structural Equation Models }
-\usage{summary.adjchisq(adj.obj)}{print.adjchisq(adj.obj)}
+\usage{summary.adjchisq(adj.obj)
+
+ \method{print}{adjchisq}(adj.obj)
+}
+
\arguments{
\item{adj.obj}{An adjusted Chi Squre object from one of the several methods meant to deal with structural equations models with non-normal data.}
}
Added: pkg/man/anova.adjchisq.Rd
===================================================================
--- pkg/man/anova.adjchisq.Rd (rev 0)
+++ pkg/man/anova.adjchisq.Rd 2010-04-21 06:38:21 UTC (rev 5)
@@ -0,0 +1,45 @@
+\name{anova.adjchisq}
+\alias{anova.adjchisq}
+\title{Adjusted Chi Squared Tests for Structural Equation Models}
+\description{ Adjusted Chi Squared Tests for Structural Equation Models \link{sbchisq}}
+
+\usage{anova.adjchisq(adjobj0, adjobj1)
+
+}
+
+\arguments{
+ \item{adjobj0}{An adjusted Satorra-Bentler Chi Squre object from \link{sbchisq}}
+ \item{adjobj1}{An adjusted Satorra-Bentler Chi Squre object from \link{sbchisq}}
+}
+
+\details{Compares two nested models that have been corrected for non-normality of the data}
+
+\value{Prints a matrix of the test similar to \link{anova}. Note, due to scaling, the Chi Square evaluated will not be the same as the difference between the two model Chi Square values.
+
+The corrected chi square difference test is calculated as follows - As one model is nested within the other, let's call m0 the nested model and m1 the more general model. DF = df0 - df1.
+
+To calculate the chi square difference between these two models, we need to look at the correction factors - c0 and c1 - and calculate a correction factor for the test itself, cd.
+
+\deqn{cd = (df0* chisq0 - df1 * chisq1)/(df0 - df1)}
+
+where chisq is the uncorrected chi square value for each model. The Chi Square statistic, TRd, is then defined as follows:
+
+\deqn{TRd = (chisq0 - chisq1)/cd}
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/sem-additions -r 5
More information about the Sem-additions-commits
mailing list