Update simple slope analysis

Posted on 2009-10-24 by

In the article ,I made basic introduction about simple slope, and the example included involve using of the provided dataset in R, from the package of QuantPsych. After that, I developed this topic to an article. In my own paper (now is under review), I substitute the data provided by R for applied research data. The trick is simulation, by applying efficient statistics (mean, sd, the correlation matrix of different variables, the most important, the value of correlaitions between interaction and other variables) provided in published academic articles. I simulated data set with almost the same relationships between variables, for anyone interested in this topic, you can try this method following my codes (in R software and ). By processing these programs, you can apply raw information contained in any research to proceed your own analysis concerning simple slopes to explore the simple effect you interested in but the author didn't analyze in their articles. Note that my paper concern cases with continuous variables as moderator and mainly focus on linear regression.
the following contain both English and Chinese version of my code
updata: the lastest codes are not available temperorily.
English version:
data simulation:click here
simple slope analysis:click here


Chinese version:
数据模拟: 点击 这里

简单斜率分析：点击这里

Posted in R, psychology, statistics | Tagged , , , | Leave a comment

Free dataset for statistic methods

Posted on 2008-12-30 by

All the work I do today(the day next to last of 200 8) is looking for free datasets for a small paper I'm working on. It's not easy because nothing comes free.However, everything comes to him who waits and finally there is something good that will happen right around the corner. Following is my sharing:
Free datasets from UCLA
From online math center
These two are most useful and there are too many links on them!

Posted in other readings, statistics | Tagged , | Leave a comment

simple slope test

Posted on 2008-12-08 by

Notice:If you want to apply the provided dataset in R to practise simple slope test, please see the latest version

simple slope test，一种检验调节效应的方法。
传统方法直接看回归方程的交互效应项的回归系数判定是否存在调节效应，如果存在交互作用且回归系数显著，则存在调节变量（属于传统的假设检验体系）。之后Johnson–Neyman，提出显著域的方法，与区间估计一致，先计算出显著水平的调节变量z的取值，利用这个值通过t统计量构造一个显著域，即在显著域（z值的某个范围）内，交互作用的回归系数是显著的。而simple slope test则在发现显著交互效应之后利用简单回归方程对调节变量如何对主要自变量与因变量的关系产生影响进行了进一步探索。所谓简单回归方程，即在调节变量z的某几个固定水平下，y对主要自变量x的回归。simple slopes就是指的在这种简单回归方程中x变量的回归系数，即不同水平调节变量x与y之间的关系。到现在，simple slope test所报告的已经不仅仅是呈现简单方程的交互作用图，以及简单方程回归系数（也就是simple slopes ）的t值，p值显著性，回归系数的区间估计已经是一个更为合理的选择。
说不如做，经过一天的探索，利用R中QuantPsyc包中的sim.slopes函数的代码经过修改，进一步延展了原代码的不足，将simple slopes的区间估计，简单方程交互作用图，以及simple slopes的significant region都包含在其中。下图即利用QuantPsyc包中提供的实例数据所做的一个simple slope test:simple slope test
图左边即简单回归方程所体现出来的交互作用，即调节作用，右图显示的调节作用的significance region之前的代码中significant region的计算错误，可对比看代码和图看错误的地方，simple slope(update)在$latex z_{\frac{\alpha}{2}$及$latex z_{1-\frac{\alpha}{2}}$之间就是没有调节效应区域，z值的其他水平就是存在调节效应发生影响的水平范围。这样看，传统回归的检验方法仅仅解释了调节效应是什么的问题主要关注调节效应存在与否极其效应大小(如果计算效应量)，而simple slope test 则进一步解释了调节效应如何产生影响的问题。R代码
见更新R代码
正确计算显著域为：
$latex Z>-\frac{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
cov(B_{1},B_{3})-B_{1}\times
B_{3}}{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{3})-B_{3}^2}+\sqrt{(\frac{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
cov(B_{1},B_{3})-B_{1}\times
B_{3}}{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{3}-B_{3}^2)})^2+\frac{B_{1}^2-t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{1})}{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{3})-B_{3}^2}}$，
或者$latex Z<-\frac{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
cov(B_{1},B_{3})-B_{1}\times
B_{3}}{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{3})-B_{3}^2}-\sqrt{(\frac{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
cov(B_{1},B_{3})-B_{1}\times
B_{3}}{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{3})-B_{3}^2})^2+\frac{B_{1}^2-t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{1})}{t_{1-\frac{\alpha}{2}(n-k-1)}^2\times
var(B_{3})-B_{3}^2}}$

Posted in R, psychology, regression | Tagged , | 6 Comments

R中计算SD时缺失变量的影响

Posted on 2008-12-06 by

R中用于进行simple slope test的 包不能够正常运用，于是打算仿照其程序自己写下，过程中发现当数据中存在缺失值时直接利用sd命令计算SD存在一定问题。运行如下命令
##a=c(runif(10),NA)
de1=sd(a,na.rm = TRUE)
de2=sd(a,na.rm = FALSE)
得到结果
> de1
[1] 0.2764684
> de2
[1] NA

na.rm命令用于决定是否删除缺失值，当选择TRUE时，得到的是以N-1计算的SD，而选择FALSE时理论上应该包含缺失值进行计算，应该得到一个以N计算的SD，但结果是返回NA。查了小半天，发现原因所在，sd函数中也是通过sqrt(var())命令来进行SD的计算 ，而在var()命令中，同样需要进行na.rm的设置，当na.rm设置为FALSE而不对 use 参数进行设置时，use 就缺省地使用"everything"这一选项进行，在这种情况下，

If use is "everything", NAs will propagate conceptually, i.e., a resulting value will be NA whenever one of its contributing observations is NA

问题就出在这里，当直接使用sd命令而数据又含有缺失值时，只能够获得删除缺失值的结果，而如果要计算包含缺失值的SD，需要间接地采用sqrt(var(,use="complete.obs"))来进行计算

Posted in R | Tagged , | 2 Comments

搞清楚的问题和不清楚的问题

Posted on 2008-05-27 by

接着上次搞清楚的GMM模型中Class变量百分比的设定，今天比较轻松地搞懂了论文（Deciding on the Number of Classes in Latent Class Analysis and Growth Mixture Modeling: A Monte Carlo Simulation Study, STRUCTURAL EQUATION MODELING, 14(4), 535–569）中各四个LCA模型的设定，并且实际生成了数据。当然在LCA中，Class 变量百分比的设定完全无差，发现论文中的一个小错误，在两个模型的百分比中作者将设定的各类别的百分比之和居然大于100%，这应该是不正确的吧。将文中模型设定部分的item probability 理解为continuous 变量的均值，而item means理解为categorical 变量的均值，如我的理解进行模型设定，生成数据的结果显示，continuous 变量的均值population确实就是item probability，而categorical 变量同样。不过一个问题是，文中的模型假定部分，item probability 对应的括号中是categorical 变量的均值，对categorical 变量来说，那个均值是threshold。我认为对于continuous变量来说均只代表的item归属于某一类别的可能性，那么对categorical 变量来说，我能够理解的是这个threshold 就像坐标轴上的某个点，但如何由 threshold 得到其归属某一类别 probability？中间的转换关系是怎样？

Posted in Mplus | 1 Comment

随机生成Class 变量的百分比设定

Posted on 2008-05-24 by

经过小半天的实验，终于搞懂了如何在class变量中设定个类别的百分比。虽然有点后知后觉，不过还是经过自己动手做才能彻底明白其中的关系。开始一直在尝试改变一个简单模型中[c#1*?]的设定值，然后看结果中latent proportion 的变化，结果发现这项结果几乎没有变化，当然也没有很清楚对于名义变量和类别变量均值的含义，因此直接将期均值定位自己想要的百分比，尝试自然以失败告终。然后直接从模拟得到的数据结果出发来看规律，一个尝试的结果：

[C#*?] Continue reading →

Posted in Mplus | 1 Comment

一例简单的spss双因素方差分析实践

Posted on 2008-05-21 by

昨天被同学问到一个问题，同学问题很简单，就是用一个summary 的数据做一个双因素方差分析，然后作出交互做用图（她的目的MS只是为了用图），数据如下：

y1 y2
x1 4 8
x2 2 3

一开始有点懵，往常经验都是用原始数据来做ANOVA，summary 数据怎么看都觉得有点奇怪。开始直接把四格表数据在Excel里边套用双因素方差分析，出来的结果很诡异，而且觉得这种逻辑相当不对，直接套数据，不清楚因变量自变量。然后好好理了思路，即使要做方差分析，不还原原始数据，四格表也应该重新进行数据表达，如果用summary 数据，下面这个才是我们应该用来进行方差分析的数据表：

表2
x1 x2 y
1 1 4
1 2 8
2 1 2
2 2 3

但利用上面数据在Excel 中进行重复/无重复的方差分析均为得到正确结果，看来在Excel里边这样的数据形式是不对的。导出在SPSS里边进行Univariate分析，利用上面的表2 summary 数据形式，当选择fix factor时不可能进行分析，因为其残差部分为0，因此不可能计算出相应F值和p-value；但选择random factor时则可以进行分析，能够得到分析结果，结果显示简单主效应均不显著，而交互作用项残差为0，依然无法进行显著性检验。
random summary
不知道问题出在哪里，于是用R产生了两列四个样本量为10，均值分别对应summary 数据的正态随机数，进行分析。当采用fix factor 时，分析结果简单主效应和交互效应均是显著的，fix factor
而采用random factor 策略对原始数据进行分析，得到以下结果：random factor
利用原始数据可以进行交互效应的分析，但在这种分析之下，简单主效应都变得不显著而只有交互效应是显著的。
从上面两种不同方法，可以直观地看出random 和fix factor 方法的区别。在random 中，利用计算F值的MSE是变量的残差，即组内残差，
“A” F = Mean SquareA / Mean SquareAxB
“B” F = Mean SquareB / Mean SquareAxB
“AxB” F = Mean SquareAxB / Mean SquareERROR

而在fixed factor方法中，利用计算F值的MSE是总体的残差，
“A” F = Mean SquareA / Mean SquareERROR
“B” F = Mean SquareB / Mean SquareERROR
“A B” F = Mean SquareAxB / Mean SquareERROR
更进一步说，在fixed factor 条件之下，作为自变量的x是看做一个恒定的量，因此其本身不存在误差，在这种模型中误差全部合并到一个误差变量中，考察变量的作用；而在random factor 中，进行分析的自变量x是看做一个随机变量，因此本身就具有一个误差项，这个误差项是代表其作为随机变量的抽样误差。其实这个里边不太明白的是为何分母是交互作用项？同时从数据来看的话，交互作用项与X变量的误差项是相等的，这个中间的关系不太清楚，有待搞清楚。

Posted in statistics | 2 Comments

简单模拟

Posted on 2008-05-19 by

以手册上的简单例子为模型，进行了简单的真模型模拟，发现以下结果：
title:
this is an example of a GMM for a
continuous outcome using automatic
starting values and random starts
montecarlo:
names are y1-y4 x;
genclasses = c(2);
classes = c(2);
nobs = 300;
seed = 3454367;
nrep = 100;
save =monte.dat;
analysis:
type = mixture;
lrtbootstrap = 100;
model population:
%overall%
[x@0]; x@1;
y1-y4*.5;
i s | y1@0 y2@1 y3@2 y4@3;
i*1; s*.2;
c#1 on x*1;
i on x*.5;
s on x*.3;
%c#1%
[i*1 s*.5];
%c#2%
[i*3 s*1];
model:
%overall%
y1-y4*.5;
i s | y1@0 y2@1 y3@2 y4@3;
i*1; s*.2;
c#1 on x*1;
i on x*.5;
s on x*.3;
%c#1%
[i*1 s*.5];
%c#2%
[i*3 s*1];
output:
tech13 tech9 tech11 tech14;
1） AIC, BIC, SABIC,在同一个模型中$latex chi^2$检验的结果中observed proportion是一样的？变换样本同样的情况？这该如何解释。我的理解是：AIC, BIC, SABIC虽然都服从$latex \chi^{2}$分布，但根据他们呢的计算公式，不可能是服从同一个分布，应该是有各自的分布，但为何结果中proportion 会是一样？（proportion即在replication中观测到的值超过在理论分布中对应百分比的那个数值的实际百分比）
Nobserve=300, Nreps=100,
AIC BIC SABIC
0.99 0.99 0.99
0.99 0.99 0.99
0.92 0.92 0.92
0.87 0.87 0.87
0.81 0.81 0.81
0.75 0.75 0.75
0.51 0.51 0.51
0.32 0.32 0.32
0.21 0.21 0.21
0.09 0.09 0.09
0.04 0.04 0.04
0.01 0.01 0.01
0.01 0.01 0.01
2）

Nobserve LMR MST MKT BLRT
100 30.00% 0.00% 0.00% 20.00%
300 50.00% 20.00% 20.00% 50.00%
500 100.00% 10.00% 0.00% 100.00%
1000 100.00% 10.00% 0.00% 100.00%
在指数部分，对于LMR 以及BLRT来说，样本量对其影响似乎有某个上限水平，例如在样本量为100的时候，LMR与BLRT正确识别的概率均为较低，当样本量达到300时二者马上上升到50%与其说受样本量的影响大不如说受重复次数的影响更大，取500样本量，时，在5%水平，他们识别出正确模型的水平都达到100%，此时样本量再上升对他们的已经不再有影响。（不过这个情况有可能因为重复次数太少，在这几个模拟中都只有10次重复）。于是在尝试了提高重复次数。
3）在样本量为300 的条件下进行了一个对比：发现在同等条件下，重复次数会影响指数的表现。但对不同指数影响不同。这一影响应该同样本量有紧密关系，但还没尝试。

NREPS LMR MST MKT BLRT
10 50.00% 20.00% 20.00% 50.00%
100 54.00% 7.00% 8.00% 63.00%

Posted in Mplus | 1 Comment

Mplus 一个实例疑问？

Posted on 2008-05-07 by

利用结果方程课程p17页的一个简单例子，先根据已知的factor loading 生成模拟数据，再利用数据拟合模型，结果发现奇怪的结果。首先根据Mplus user guide p53ex5.1中的例子的Monte Carlo 模拟例子生成数据。

title: this is an example of a CFA model
montecarlo:
names = y1-y9;
nobs = 500;
nreps = 1;
save = sem1.dat;
model montecarlo:
f1 by y1@.73 y4@.69 y5@.65 y2@.17 y6@.16;
f2 by y7@.14 y3@.66 y8@.80 y9@.67;
f1 with f2*.24;
f1-f2*1;
y1-y9*1;
model:
f1 by y1@.73 y4@.69 y5@.65 y2@.17 y6@.16;
f2 by y7@.14 y3@.66 y8@.80 y9@.67;
f1 with f2*.24;
y1-y9*1;
f1-f2*1;
output:
tech9;
然后根据CFA的语法程序来拟合数据。
TITLE: this is an example of a CFA model
DATA: FILE IS sem1.dat;
VARIABLE: NAMES ARE y1-y9;
MODEL: f1 BY y1 y4 y5 y2 y6;
f2 BY y7 y3 y8 y9;
结果很奇怪：
MODEL RESULTS
Estimates
F1 BY
Y1 1.000
Y4 0.915
Y5 1.008
Y2 0.317
Y6 0.271
F2 BY
Y7 1.000
Y3 26.562
Y8 33.555
Y9 29.181
F2 WITH
F1 0.000
Variances
F1 0.510
F2 0.000
Residual Variances
Y1 0.982
Y2 1.033
Y3 1.129
Y4 1.065
Y5 1.015
Y6 0.985
Y7 0.999
Y8 0.801
Y9 1.005
除了f1 的loading 部分和Y的残差部分还算正常外，其他的拟合差距都很大；
而如果生成数据的时候将f1-f2的方差设定为0，进行数据的生成，再用同样的模型来拟合，效果就很改进很多。
MODEL RESULTS
Estimates S.E. Est./S.E.
F1 BY
Y1 1.000 0.000 0.000
Y4 1.142 0.636 1.795
Y5 1.854 0.879 2.108
Y2 1.332 0.698 1.910
Y6 0.909 0.562 1.617
F2 BY
Y7 1.000 0.000 0.000
Y3 2.502 1.269 1.972
Y8 1.565 0.832 1.881
Y9 1.814 0.953 1.904
F2 WITH
F1 0.027 0.016 1.657
Variances
F1 0.008 0.011 0.712
F2 -0.012 0.012 -1.075
Residual Variances
Y1 0.965 0.062 15.675
Y2 1.027 0.067 15.402
Y3 1.104 0.088 12.566
Y4 1.021 0.065 15.595
Y5 1.033 0.072 14.373
Y6 0.984 0.063 15.727
Y7 1.002 0.065 15.371
Y8 0.868 0.060 14.429
Y9 1.029 0.072 14.210
在手册的例子中，在生成数据和进行模型拟合时用完全相同的Model，结果也不差。所以问题是：真的方差的缘故吗？
如果更复杂一点的SEM模型，例如3个factor,在生成数据的Model部分，通常会出现POPULATION COVARIANCE的错误，例如，假设模型3个因子间存在关系f1 with f2， f2 with f3 ，f1 with f3，在生成数据时，如果3个关系保留，就会出现上述错误，而如果删除f1 with f3的关系，就正确；删除其他任何一个关系错误依然存在；所以到现在也不理解问题的关键？
另外小的tips 在于，在某些场合下，Mplus对语句出现的先后顺序有严格要求；对于语法中= () 与变量之间必须有空格；某个语句结束时候一定要;……总的说来,Mplus的程序语法似乎比lisrel要简洁些，不过某些语法还是比较诡异。

Posted in Mplus | 8 Comments

To understand CI

Posted on 2007-12-11 by

Here is a R-script for simulating the process of computing CI for given samples
n: the mumber of samples
###:sz: the sample size of every sample
###:S:the sample
###:alpha:the significant level
f = function(n, alpha, sz ) {
S= matrix( rnorm(sz*n),ncol=n,nrow=sz)
m =c( colMeans(S))
z = qnorm(1 - alpha/2)
y0 = m - z * 1/sqrt(sz)
y1 = m + z * 1/sqrt(sz)
plot(1, xlim = c(0.5, n + 0.5), ylim = c(min(y0), max(y1)),
type = "n", xlab = "Samples", ylab = "Confidence interval",
main = expression("CI: [" ~ mu - z[alpha/2] * sigma/sqrt(n) ~
", " ~ mu + z[alpha/2] * sigma/sqrt(n) ~ "]"))
abline(h = 0, lty = 2)
for (i in 1:n) {
arrows(i, y0[i], i, y1[i], length = 0.05, angle = 90,
code = 3, col = ifelse(0 > y0[i] & 0 < y1[i], "gray", "red")) points(i, m[i]) } } f(50, 0.05, 100)

From the graph we can try to understand the true meaning of confidence level, "The probability that we are correct.It was the frequency of correct statments that a statistician sho use his method will make in the long run.It says noting about how'accurate the current estimate is."(《Lady Tasting Tea》,p123),the red ones in the graph is the 5% don't include $latex \mu$
PS: A question to Mr.Li: in the paper Inference by Eye:Confidence Intervals and How to Read Pictures of Data (Geoff Cumming,Sue Finch,2005) in Pge 175,it's said that"A large w may suggest an experiment is of little value", I didn't makemuch sense of this sentence. I'm not sure if it's right to understand it this way: because $latex w=t_{(n-1),c}\times\frac{SD}{\sqrt{n}}$,given n and $latex \alpha$, only SD could exerts effect on w, so there must be a large SD, morever, large SD imlpies errors in large scale,so large w indicate bad experimental control, and finally resoned the conclusion that "is of little value"

--
David Salsburg(2001), The Lady Tasting Tea:how statistics revolutionized science in the twentieth century,W.H.Freeman AND Company New York

http://www.yihui.name/en/index.php

Posted in R | Leave a comment