《结构方程模型及其应用》(侯, 温, & 成,2004)部分章节R代码

John FOX教授的sem包试写了这本教材的几个例子,结果都与LISREL8报告的Minimum Fit Function Chi-Square吻合。不过LISREL其它拟合指标用的是Normal Theory Weighted Least Squares Chi-Square?,所以看上去比Minimum Fit Function Chi-Square报告的结果要好那么一些些。LISREL历史上推的GFI/AGFI曾因为经常将差报好被批评,圈内朋友私下嘲笑这样LISREL就更好卖了,买软件的用户也高兴将差报好,只有读Paper的人上当。

目前只写了chap3_1..到Chap3_2_...,一共5(或6)个例子。尺有所短,寸有所长--sem包不会自动报告所有修正指数,不能做样本量不一样多的多组模型,对复杂的模型要写的代码太多。不过,在已经尝试的几个例子里,有一个是LISREL跑不出来但sem包能跑出结果的。目前sem包还没有达到结构方程众多商业软件的成熟水准,但R庞大的义工武器库已经使sem包至少已经在Missing Data Multiple Imputation、Bootstrapping等等应用上胜出一筹。所以例子中还附带了缺失数据一讲Multiple Impuation的R示范代码。

欢迎各位有兴趣作类似尝试的同学将结果email我, 可以陆续更新到下面的GPL版权代码集合中。

代码下载:lixiaoxu.googlepages.com (中大镜像)

[update] AMos, Mplus 与 SAS/STAT CALIS 缺省报告与使用的都是 Minimum Fit Function Chi-Square,可通过各种软件(Albright & Park, 2009)的结果对比查验。Normal Theory Weighted Least Squares Chi-Square并非总是比Minimum Fit Function Chi-Square报告更“好”的拟合结果(虽然常常如此)。Olsson, Foss, 和 Breivik (2004) 用模拟数据对比了二者,确证Minimum Fit Function Chi-Square计算得到的拟合指标在小样本之外的情形都比Normal Theory Weighted Least Squares Chi-Square的指标更适用。

R的sem包目前在迭代初值的计算上还做得不够好。我遇到的几个缺省初值下迭代不收敛案例,将参数设定合理范围的任意初值(比如TX,TY都设成0-1之间的正实数)之后都收敛了。

--

Albright, J. J., & Park., H. M. (2009). Confirmatory factor analysis using Amos, LISREL, Mplus, and SAS/STAT CALIS. The University Information Technology Services Center for Statistical and Mathematical Computing, Indiana University. Retrieved July 7, 2009, from http://www.indiana.edu/~statmath/stat/all/cfa/cfa.pdf.

Olsson, U. H., Foss, T., & Breivik, E. (2004). Two equivalent discrepancy functions for maximum likelihood estimation: Do their test statistics follow a non-central chi-square distribution under model misspecification? Sociological Methods Research, 32(4), 453-500.

[update]R中的sem包妙处之一是可在线实现结构方程应用界面。下面这个最粗陋的例子对应原书Chap3_1_4_CFA_MB.LS8的结果。
--荣耀属于sem package的作者Rweb的作者、以及服务器运算资源的提供者。

##Input Correlation Matrix
R.DHP<-matrix(0,ncol=17,nrow=17);
R.DHP[col(R.DHP) >= row(R.DHP)] <- c(
1,
.34,1,
.38,.35,1,
.02,.03,.04,1,
.15,.19,.14,.02,1,
.17,.15,.20,.01,.42,1,
.20,.13,.12,.00,.40,.21,1,
.32,.32,.21,.03,.10,.10,.07,1,
.10,.17,.12,.02,.15,.18,.23,.13,1,
.14,.16,.15,.03,.14,.19,.18,.18,.37,1,
.14,.15,.19,.01,.18,.30,.13,.08,.38,.38,1,
.18,.16,.24,.02,.14,.21,.21,.22,.06,.23,.18,1,
.19,.20,.15,.01,.14,.24,.09,.24,.15,.21,.21,.45,1,
.18,.21,.18,.03,.25,.18,.18,.18,.22,.12,.24,.28,.35,1,
.08,.18,.16,.01,.22,.20,.22,.12,.12,.16,.21,.25,.20,.26,1,
.12,.16,.25,.02,.15,.12,.20,.14,.17,.20,.14,.20,.15,.20,.50,1,
.20,.16,.18,.04,.25,.14,.21,.17,.21,.21,.23,.15,.21,.22,.29,.41,1
);
R.DHP<-t(R.DHP);
colnames(R.DHP)<-rownames(R.DHP)<-paste('X',1:17,sep='');
print('Inputted Correlation Matrix');
print(R.DHP);
##
##
##Input Model Specification of Chap3_1_4_CFA_MB.LS8
require(sem);
model.B <- matrix(ncol=3,byrow=TRUE,data=c(
'X1 <-> X1' , 'TD1_1'  , NA  ,
'X2 <-> X2' , 'TD2_2'  , NA  ,
'X3 <-> X3' , 'TD3_3'  , NA  ,
'X5 <-> X5' , 'TD5_5'  , NA  ,
'X6 <-> X6' , 'TD6_6'  , NA  ,
'X7 <-> X7' , 'TD7_7'  , NA  ,
'X8 <-> X8' , 'TD8_8'  , NA  ,
'X9 <-> X9' , 'TD9_9'  , NA  ,
'X10<-> X10', 'TD10_10', NA  ,
'X11<-> X11', 'TD11_11', NA  ,
'X12<-> X12', 'TD12_12', NA  ,
'X13<-> X13', 'TD13_13', NA  ,
'X14<-> X14', 'TD14_14', NA  ,
'X15<-> X15', 'TD15_15', NA  ,
'X16<-> X16', 'TD16_16', NA  ,
'X17<-> X17', 'TD17_17', NA  ,
'xi1<-> xi1', NA       , '1' ,
'xi2<-> xi2', NA       , '1' ,
'xi3<-> xi3', NA       , '1' ,
'xi4<-> xi4', NA       , '1' ,
'xi5<-> xi5', NA       , '1' ,
'xi1<-> xi2', 'PH12'   , NA  ,
'xi1<-> xi3', 'PH13'   , NA  ,
'xi1<-> xi4', 'PH14'   , NA  ,
'xi1<-> xi5', 'PH15'   , NA  ,
'xi2<-> xi3', 'PH23'   , NA  ,
'xi2<-> xi4', 'PH24'   , NA  ,
'xi2<-> xi5', 'PH25'   , NA  ,
'xi3<-> xi4', 'PH34'   , NA  ,
'xi3<-> xi5', 'PH35'   , NA  ,
'xi4<-> xi5', 'PH45'   , NA  ,
'X1 <- xi1' , 'LX1_1'  , NA  ,
'X2 <- xi1' , 'LX2_1'  , NA  ,
'X3 <- xi1' , 'LX3_1'  , NA  ,
'X5 <- xi2' , 'LX5_2'  , NA  ,
'X6 <- xi2' , 'LX6_2'  , NA  ,
'X7 <- xi2' , 'LX7_2'  , NA  ,
'X8 <- xi1' , 'LX8_1'  , NA  ,
'X9 <- xi3' , 'LX9_3'  , NA  ,
'X10<- xi3' , 'LX10_3' , NA  ,
'X11<- xi3' , 'LX11_3' , NA  ,
'X12<- xi4' , 'LX12_4' , NA  ,
'X13<- xi4' , 'LX13_4' , NA  ,
'X14<- xi4' , 'LX14_4' , NA  ,
'X15<- xi5' , 'LX15_5' , NA  ,
'X16<- xi5' , 'LX16_5' , NA  ,
'X17<- xi5' , 'LX17_5' , NA  )
);
class(model.B)<-'mod';
##
##
N=350;##sample size;
##R.DHP[-4,-4] excludes X4
## Result
(summary(sem.B<-sem(model.B, R.DHP[-4,-4], N)));
## Residuals
(round(residuals(sem.B),3));
#####################
boxplot.matrix = function(M,ylim=c(-1,1)) {
  M = as.matrix(M);
  boxplot(c(M[row(M)>col(M)]),at=1,xlab='',ylab='',ylim=ylim);
  points(rep(1,length(c(M[row(M)>col(M)]))),c(M[row(M)>col(M)]),pch='-',col='red');
  stem(c(M[row(M)>col(M)]));
  boxplot.stats(c(M[row(M)>col(M)]));
}
####################
boxplot(residuals(sem.B));

Paper for 1st Chinese useR! Conference: Web Powered by R, or R Powered by Web

欢迎在本部的同学明天上午到现场看李崇亮同学演示,地点见会议主页

论文下载(Googlepages, 中文大学镜像)
RWebFriend for WordPress 在线示例(yo2.cn上的示例,  奇想录上的临时示例)

Google Presentation 在线演示

[update2009.07.11]文中MediaWiki与R集成的实验平台已不再开放编辑权限,原例转到另一个平台上(界面是西班牙语,找不到英语界面的平台)。希望能有帮手提供linux服务器自建一个平台。

“Confidence interval of R-square”, but, which one?

In linear regression, confidence interval (CI) of population DV is narrower than that of predicted DV. With the assumption of generalizability, CI of at is

,

while CI of is

.

The pivot methods of both are quite similar as following.

,

so .

,

so

of linear regression is the point estimate of

for fixed IV(s) model. Or, it is the point estimate of wherein denotes the correlation of Y and , the linear composition of random IV(s) . The CI of is wider than that of with the same and confidence level.

[update] It is obvious that CI of  relies on the distribution presumption of IV(s) and DV, as fixed IV(s) are just special cases of generally random IV(s). Usually, the presumption is that all IV(s) and DV are from multivariate normal distribution.

In the bivariate normal case with a single random IV, through Fisher's z-transform of Pearson's r, CI of the re-sampled can also be constructed. Intuitively, it should be wider than CI of .

Thus,

CI of can be constructed as . With the reverse transform , the CI bounds of are

and

.

In multiple p IV(s) case, Fisher's z-transform is

.

Although it could also be used to construct CI of , it is inferior to noncentral F approximation of R (Lee, 1971). The latter is the algorithm adopted by MSDOS software R2 (Steiger & Fouladi, 1992) and R-function ci.R2(...) within package MBESS (Kelley, 2008).

In literature, "CI(s) of R-square" are hardly the literal CI(s) of in replication once more. Most of them actually refer to CI of . Authors in social science unfamiliar to hate to type when they feel convenient to type r or R. Users of experimentally designed fixed IV(s) should have reported CI of . However, if they were too familiar to Steiger's software R2 to ignore his series papers on CI of effect size, it would be significant chance for them to report a loose CI of , even in a looser name "CI of ".

----

Lee, Y. S. (1971). Some results on the sampling distribution of the multiple correlation coefficient. Journal of the Royal Statistical Society, B, 33, 117–130.

Kelley, K. (2008). MBESS: Methods for the Behavioral, Educational, and Social Sciences. R package version 1.0.1. [Computer software]. Available from http://www.indiana.edu/~kenkel

Steiger, J. H., & Fouladi, R. T. (1992). R2: A computer program for interval estimation, power calculation, and hypothesis testing for the squared multiple correlation. Behavior research methods, instruments and computers, 4, 581–582.

R Code of Part I:



R Code of Part II:



WordPress (and WPMU) Plugin for R Web Interface

Download: RwebFriend.zip [Update] Including Chinese UTF8 Version
Plugin Name: RwebFriend
Plugin URL: http://xiaoxu.lxxm.com/RwebFriend
Description: Set Rweb url options and transform [rcode]...[/rcode] or <rcode>...</rcode> tag-pair into TEXTAREA which supports direct submit to web interface of R.

*Credit notes:codes of two relevant plugins are studied and imported. One of the plugins deals with auto html tags within TEXTAREA tag-pair, the other stops WordPress to auto-transform quotation marks.

Version: 1.0
Author: Xiaoxu LI
Author URI: http://xiaoxu.lxxm.com/
Setup:Wordpress 3.5

ScreenShot

 

WordPress 3.4

Usage:

[update] The free Chinese wordpress platform yo2.cn has installed this plugin.  See my demo.

More online demos -- http://wiki.qixianglu.cn/rwebfriendttest/

[update, June 2009] 72pines.com (here!) installed this plugin. Try



[update2009JUL18]Test installed packages of Rweb:
http://pbil.univ-lyon1.fr/cgi-bin/Rweb/Rweb.cgi
https://rweb.stat.umn.edu/cgi-bin/Rweb/Rweb.cgi

Type III ANOVA in R

Type III ANOVA SS for factor A within interaction of factor B is defined as , wherein A:B  is the pure interaction effect orthogonal to main effects of A, B, and intercept. There are some details in R to get pure interaction dummy IV(s).

Data is from SAS example PROC GLM, Example 30.3: Unbalanced ANOVA for Two-Way Design with Interaction

##
##Data from http://www.otago.ac.nz/sas/stat/chap30/sect52.htm
##
drug <- as.factor(c(t(t(rep(1,3)))%*%t(1:4))); ##Factor A
disease <- as.factor(c(t(t(1:3)) %*% t(rep(1,4))));##Factor B
y <- t(matrix(c(
42 ,44 ,36 ,13 ,19 ,22
,33 ,NA ,26 ,NA ,33 ,21
,31 ,-3 ,NA ,25 ,25 ,24
,28 ,NA ,23 ,34 ,42 ,13
,NA ,34 ,33 ,31 ,NA ,36
,3 ,26 ,28 ,32 ,4 ,16
,NA ,NA ,1 ,29 ,NA ,19
,NA ,11 ,9 ,7 ,1 ,-6
,21 ,1 ,NA ,9 ,3 ,NA
,24 ,NA ,9 ,22 ,-2 ,15
,27 ,12 ,12 ,-5 ,16 ,15
,22 ,7 ,25 ,5 ,12 ,NA
),nrow=6));
## verify data with http://www.otago.ac.nz/sas/stat/chap30/sect52.htm
(cbind(drug,disease,y));
##
## make a big table
y <- c(y);
drug <- rep(drug,6);
disease <- rep(disease,6);
##
## Design the PURE interaction dummy variables
m <- model.matrix(lm(rep(0,length(disease)) ~ disease + drug +disease:drug));
##! If lm(y~ ...) is used, the is.na(y) rows will be dropped. The residuals will be orthogonal to observed A, & B rather than designed cell A & B. It will be Type II SS rather than Type III SS.
c <- attr(m,"assign")==3;
(IV_Interaction <-residuals( lm(m[,c] ~ m[,!c])));
##
## verify data through type I & II ANOVA to http://www.otago.ac.nz/sas/stat/chap30/sect52.htm
## Type I ANOVA of A, defined by SS_A --
anova(lm(y~drug*disease));
##
## Type II ANOVA of A, defined by SS_{A+B}-SS_B --
require(car);
Anova(lm(y~drug*disease),type='II');
anova(lm(y~disease),lm(y~drug + disease))
##
##
## Type III ANOVA of A defined by SS_{A:B+A+B}-SS_{A:B+B}
t(t(c( anova(lm(y~IV_Interaction+disease),lm(y~disease * drug))$'Sum of Sq'[2]
,anova(lm(y~IV_Interaction+drug),lm(y~disease*drug))$'Sum of Sq'[2]
,anova(lm(y~disease+drug),lm(y~disease*drug))$'Sum of Sq'[2])))
##
##

Currently, Anova(...) of Prof John Fox's car package (V. 1.2-8 or 1.2-9) used "impure" interaction dummy IV(s), which made its type III result relying upon the order of factor levels. I think in its next version, the "pure" interaction dummy IV(s) will be adopted to give consistent type III SS.

[update:]

In Prof John FOX's car package, with parameter contrasts in inputted lm object, Example(Anova) gave type III SS consistent to  other softwares. In this case, the code line should be --

Anova(lm(y~drug*disease, contrasts=list(drug=contr.sum, disease=contr.sum)),type='III');

Contrasts patterns are defined within lm(...) rather than Anova(...). An lm object with default contrasts parameter is inappropriate to calculate type III SS, or the result will rely on the level names in any nominal factor --

require(car);
M2<-Moore;
M2$f1<-M2$fcategory;
M2$f2<-as.factor(- as.integer(M2$fcategory));
mod1<-lm(formula = conformity ~ f1 * partner.status,data=M2);
mod2<-lm(formula = conformity ~ f2 * partner.status,data=M2);
c(Anova(mod1,type='III')$'Sum Sq'[3],Anova(mod2,type='III')$'Sum Sq'[3])

There was hot discussion of type III ANOVA on R-help newsgroup. Thomas Lumley thought Types of SS nowadays don't have to make any real sense --

http://tolstoy.newcastle.edu.au/R/help/05/04/3009.html

This is one of many examples of an attempt to provide a mathematical answer to something that isn't a mathematical question.

As people have already pointed out, in any practical testing situation you have two models you want to compare. If you are working in an interactive statistical environment, or even in a modern batch-mode system, you can fit the two models and compare them. If you want to compare two other models, you can fit them and compare them.

However, in the Bad Old Days this was inconvenient (or so I'm told). If you had half a dozen tests, and one of the models was the same in each test, it was a substantial saving of time and effort to fit this model just once.

This led to a system where you specify a model and a set of tests: eg I'm going to fit y~a+b+c+d and I want to test (some of) y~a vs y~a+b, y~a+b vs y~a+b+c and so on. Or, I want to test (some of) y~a+b+c vs y~a+b+c+d, y~a+b+d vs y~a+b+c+d and so on. This gives the "Types" of sums of squares, which are ways of specifying sets of tests. You could pick the "Type" so that the total number of linear models you had to fit was minimized. As these are merely a computational optimization, they don't have to make any real sense. Unfortunately, as with many optimizations, they have gained a life of their own.

The "Type III" sums of squares are the same regardless of order, but this is a bad property, not a good one. The question you are asking when you test "for" a term X really does depend on what other terms are in the model, so order really does matter. However, since you can do anything just by specifying two models and comparing them, you don't actually need to worry about any of this.

-thomas

R: str(…) 与 getS3method(…,…)

感谢R专家XIE Yihui同学在线答疑

me: 请教两个R的技术:1.R中有没有对象浏览器之类的工具?一举看完一个对象的子子孙孙 2.怎么看深入的源代码> prcomp
function (x, ...)
UseMethod("prcomp")
<environment: namespace:stats>

Yihui: 1. str()是很常用的一个函数,它可以充分查看对象的子子孙孙 2. 很多函数要么是S3 method,要么是调用C code,所以一般不能直接看源代码

S3 method可以用getS3method()去查看,比如prcomp就是S3方法,那么可以看它的default方法是什么:

> getS3method('prcomp','default')
function (x, retx = TRUE, center = TRUE, scale. = FALSE, tol = NULL,
...)
{
x <- as.matrix(x)
x <- scale(x, center = center, scale = scale.)
cen <- attr(x, "scaled:center")
sc <- attr(x, "scaled:scale")
if (any(sc == 0))
stop("cannot rescale a constant/zero column to unit variance")
s <- svd(x, nu = 0)
s$d <- s$d/sqrt(max(1, nrow(x) - 1))
if (!is.null(tol)) {
rank <- sum(s$d > (s$d[1L] * tol))
if (rank < ncol(x)) {
s$v <- s$v[, 1L:rank, drop = FALSE]
s$d <- s$d[1L:rank]
}
}
dimnames(s$v) <- list(colnames(x), paste("PC", seq_len(ncol(s$v)),
sep = ""))
r <- list(sdev = s$d, rotation = s$v, center = if (is.null(cen)) FALSE else cen,
scale = if (is.null(sc)) FALSE else sc)
if (retx)
r$x <- x %*% s$v
class(r) <- "prcomp"
r
}
<environment: namespace:stats>

但所有的源代码都可以从R的源代码包中看到,Windows下经过编译了,可能有一些看不到了

me: 多谢,节约俺好多搜索时间

Yihui: 嗯,我也是花了很长时间才明白S3 method的意思,呵呵

me: 如果要看biplot.prcomp呢?

有help没有源代码的提示

Yihui: help会告诉你biplot也是generic function,可以应用在prcomp这种class上,所以:

> getS3method('biplot','prcomp')
function (x, choices = 1:2, scale = 1, pc.biplot = FALSE, ...)
{
if (length(choices) != 2)
stop("length of choices must be 2")
if (!length(scores <- x$x))
stop(gettextf("object '%s' has no scores", deparse(substitute(x))),
domain = NA)
if (is.complex(scores))
stop("biplots are not defined for complex PCA")
lam <- x$sdev[choices]
n <- NROW(scores)
lam <- lam * sqrt(n)
if (scale < 0 || scale > 1)
warning("'scale' is outside [0, 1]")
if (scale != 0)
lam <- lam^scale
else lam <- 1
if (pc.biplot)
lam <- lam/sqrt(n)
biplot.default(t(t(scores[, choices])/lam), t(t(x$rotation[,
choices]) * lam), ...)
invisible()
}
<environment: namespace:stats>

这个代码里面你会看到其实是用biplot.default,所以你要继续看default的是什么代码


Understanding QQ plots

## Try distributions like rchisq, rt, runif, rf to view its heavy, or light, left, or right tail.


n <- 30;
ry <- rnorm(n);
qqnorm(ry);qqline(ry);
max(ry)
min(ry)
##view and guess what are x(s) and y(s)
I <- rep(1,n);
qr <- ((ry%*%t(I) > I %*% t(ry))+.5*(ry %*% t(I) == I%*%t(ry)))%*%I *(1/n);##qr are the sample quantiles
points(qr,ry,col="blue"); ##to view the fact, try the following
points(qr,qr*0,col="green",pch="|");
rx <- qnorm(qr);
points(rx,ry,col="red",pch="O");
##Red O(s) circle black o(s) exactly.

03DEC2007 R-workshop sponsored by dept of psy, ZSU(=SYSU, Guang-Zhou)

Here is the updated PPT for the talk in the afternoon--which includes the zipped example codes and set-up steps for the workshop in the evening within the 3rd page. The listed anonymous on-line test (result statistics) on p-value interpretation was cited indirectly From Gigerenzer, Krauss, & Vitouch (2004).

There is an advert on http://www.psy.sysu.edu.cn/detail_news.asp?id=258 and a formal CV of the speaker is available on http://lixiaoxu.googlepageS.com

Classic Neyman-Pearson approach demo

It notes here that N-P approach does not utilize the information in the accurate p value. Actually, at the time N-P approach was firstly devised, the accurate p value was not available usually. Now almost all statistic softwares provide accurate p values and the N-P approach becomes obsolete. Wilkinson & APA TFSI (1999) recommended to report the accurate p value rather than just significance/insignificance, unless p is smaller than any meaningful precision.




--

Wilkinson, L. & APA TFSI (1999). Statistical methods in psychology journals: Guidelines and explanations. American Psychologist, 54, 594-604.