10001036 : Teaching notes of szpku dot lixiaoxu at gmail dot com

“Effect Size” — same data, different interpretations

d<-32; ## Try d<-20 ! ## to reduce the death rate by d%, From (50+d/2)% to (50-d/2)% y<-c(rep("Live",50+d/2),rep("Death",50-d/2)); y<-c(y,rep("Live",50-d/2),rep("Death",50+d/2)); y<-(y=="Live"); ## TRUE vs FALSE x<-c(rep("Treatment",100),rep("Control",100)); x<-(x=="Treatment"); ## correlation^2 cor(x,y,method="pearson")^2 cor(x,y,method="spearman")^2 cor(x,y,method="kendall")^2 ## R^2 of linear regression with norminal IV ## should we use logistic regression? ## However, R^2 is not available in GLM. ## summary(lm(y~x)) ## names(summary(lm(y~x))) summary(lm(y~x))$r.squared ## anova ## anova(lm(y~x)) ## names(anova(lm(y~x))) s<-anova(lm(y~x))$"Sum Sq";s[1]/sum(s)

Just a short R-script note to embody the 3-page-paper of Rosenthal & Rubin (1982).

Table 1. (p. 167) listed a simple set-up. There was a between-subject treatment. Control group includes 34 alive cases and 66 dead cases. Treatment group includes 66 alive cases and 34 dead cases. The question is what is the percentage of the variance explained by the nominal IV indicating the group?

The authors pointed out that one may interpret the data result as death rate was reduced by 32% while the other may interpret the same as 10.24% variance was explained. Let’s demo it more dramatically to imagine just 4% explained variance would reduce death rate by 20%.

–

Rosenthal, R. & Rubin, D. B. (1982). A simple, general purpose display of magnitude of experimental effect. Journal of Educational Psychology, 74, 166-169.

Tagged effect size, p value, R

2007 11 21 lixiaoxu English Comments (6) Permalink

Anscombe’s 4 Regressions — A Trivially Updated Demo

##---------- ## This is a trivially updated version based on the R document "?anscombe". require(stats); require(graphics) anscombe

##-- now some "magic" to do the 4 regressions in a loop:##< - ff = y ~ x for(i in 1:4) { ff[2:3] = lapply(paste(c(”y”,”x”), i, sep=”"), as.name) assign(paste(”lm.”,i,sep=”"), lmi <- lm(ff, data= anscombe)) }
## See how close they are (numerically!) sapply(objects(pattern=”lm\\.[1-4]$”), function(n) coef(get(n))) lapply(objects(pattern=”lm\\.[1-4]$”), function(n) coef(summary(get(n)))) ## Now, do what you should have done in the first place: PLOTS op <- par(mfrow=c(4,3),mar=.1+c(4,4,1,1), oma= c(0,0,2,0)) for(i in 1:4) { ff[2:3] <- lapply(paste(c(”y”,”x”), i, sep=”"), as.name) plot(ff, data =anscombe, col=”red”, pch=21, bg = “orange”, cex = 1.2, xlim=c(3,19), ylim=c(3,13)) abline(get(paste(”lm.”,i,sep=”")), col=”blue”) plot(lm(ff, data =anscombe),which=1,col=”red”, pch=21, bg = “orange”, cex = 1.2 ,sub.caption=”",caption=”" ) plot(lm(ff, data =anscombe),which=2,col=”red”, pch=21, bg = “orange”, cex = 1.2 ,sub.caption=”",caption=”" ) } mtext(”Anscombe’s 4 Regression data sets”, outer = TRUE, cex=1.5) par(op) ##
## Anscombe, F. J. (1973). Graphs in statistical analysis. American Statistician, 27, 17–21.

Tagged qq plot, R, regression

2007 11 18 lixiaoxu English Comments (1) Permalink

自由度的几何：对截距项投影残差向量的长度平方

这是《相关系数的几何：对截距投影的残差向量之间交角余弦》示意图，恰好可以用于解释为什么 $\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}$ 满足的 $\chi^2$ 分布df是n-1而不是n。

其中 $X_{i}\equiv\mu+\varepsilon_{i}$ 且 $\left[\begin{array}{c}\varepsilon_{1}\\\varepsilon_{2}\\\vdots\\\varepsilon_{n}\end{array}\right]$ 是n维空间中的标准正态随机向量。那么，容易知道有 $\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}=\sum_{i=1}^{n}\left(\varepsilon{}_{i}-\bar{\varepsilon}\right)^{2}$ 。这个表达式就是向量 $\left[\begin{array}{c}\varepsilon_{1}\\\varepsilon_{2}\\\vdots\\\varepsilon_{n}\end{array}\right]-\left[\begin{array}{c}\bar{\varepsilon}\\\bar{\varepsilon}\\\vdots\\\bar{\varepsilon}\end{array}\right]$ 长度的平方。我们已经知道， $\left[\begin{array}{c}\bar{\varepsilon}\\\bar{\varepsilon}\\\vdots\\\bar{\varepsilon}\end{array}\right]$ 就是 $\left[\begin{array}{c}\varepsilon_{1}\\\varepsilon_{2}\\\vdots\\\varepsilon_{n}\end{array}\right]$ 在截距向量（日晷指针） $\left[\begin{array}{c}1\\1\\\vdots\\1\end{array}\right]$ 上的投影。自然， $\left[\begin{array}{c}\varepsilon_{1}\\\varepsilon_{2}\\\vdots\\\varepsilon_{n}\end{array}\right]-\left[\begin{array}{c}\bar{\varepsilon}\\\bar{\varepsilon}\\\vdots\\\bar{\varepsilon}\end{array}\right]$ 就是对截距项投影残差向量，也就是在日晷盘上的投影。

日晷所处空间的n是3。如果我们对 $\left[\begin{array}{c}\varepsilon_{1}\\\varepsilon_{2}\\\varepsilon_{3}\end{array}\right]$ 抽样许多次，就会看到三维空间中各个方向对称的标准正态分布散点图。这些散点图在日晷盘上的投影就是二维空间标准正态分布散点图。日晷盘中这些点对应向量的长度平方自然是 $\chi^2_{df=2}$ 的抽样。

Tagged correlation, geometry

2007 11 16 lixiaoxu Chinese Comments (0) Permalink

习题：一类错误的注水

一个研究者每次都先看一下计算出的统计量再决定对零假设 $\mu=0$ 做单尾检验还是双尾检验。如果统计量 $\bar{X}>0$ ，就设对立假设为 $\mu>0$ ；如果统计量 $\bar{X}<0$ ，就设对立假设为 $\mu<0$ 。假如他的 $\alpha=0.05$ 请问他真实的一类错误率是多少？具体说，有许多次的实验，真实情形都是 $\mu=0$ ，他能检验出显著拒绝的比例会趋近于多少？

Tagged Hypothesis Test

2007 11 14 lixiaoxu Chinese Comments (4) Permalink

读汪丁丁《为中国股民找一个理由》所想到与读到

汪文注明首发《IT经理世界》，我读于CCER新闻。特别摘出部分：

…当代中国人生活在千年未有之变局之内，经历着三重转型同时发生的阶段，所以，每一中国人的未来，充满着奈特意义上的不确定性。这种不确定性是不可预期且不可重复的，当代实验经济学家称之为“ambiguity”，以区分于“risk(风险)”。

股票市场固然风险很高，可是，我们综观股市之外的种种生活，风险何尝不高呢？与其走出股市得一人生之幻灭，为何不走入股市搏一幻灭之人生？或者，用经济学的术语再说一遍：股市之外的高风险人生，却并不带来相应的高回报。大众纷纷进入股市，是因为他们知道在那里承担风险至少有带来相应回报的可能性。也就是说，与其终生储蓄在银行里并希望渺茫地预期不断上涨的养老、医疗、住房、教育和日常生活的费用不至于完全侵蚀了他们微不足道的储蓄，不如以这微不足道的储蓄充当投资股市的本钱，反而是更富理性的选择。

…

这涉及到我备课过程中原先没意识到也许密切关联的两个论题。第一是奈特/Knight(1921)的可测度的Risk和不可测度的Uncertainty的区分。汪文第二段中的“风险”显然是Uncertainty而不是Risk。有意思的是，不可测度的Uncertainty却是可比较高低的（这不是汪的创见，而是Knight原著的意见）。用心理计量学术语，Uncertainty不是scale变量，但却是ordinal变量，而且很可能还是连续的ordinal变量。

Knight原著并不易读。甚至只是翻查《新帕尔格雷夫经济学大辞典》1987版1996中译本的Uncertainty和Knight辞条，就已经令人云里雾里。其中Knight辞条执笔者是G. J. Stigler，他对Knight在Uncertainty上的“贡献”略有微词。Knight原著第7章注解1也小心的指出他打算规避认识论/知识论的讨论。这给我的感觉就好比：讨论一个被定义为“本质上不可讨论的对象”的对象。须知Uncertainty在Knight原著中唯一的内涵就是不可测度，于是所有对它的减少(eliminate)都是对它的否定。一旦比较它有多么地“不可测度”，就是在否定“不可测度”的本质。从罗素悖论的经验，我实在怀疑“不可测度性”程度的比较注定要引出悖论。

这便引出与之相联系的第二个论题：“主观概率”。在Uncertainty辞条中Knight的角色只是一笔带过，而主观概率才是更实质的关键词。似乎很根本的一个问题是：如果我们“完全地、本质地”不知道一个随机分布，在何种程度上能或者不能建立起一个有普遍意义的主观概率分布？–也许读懂辞条后，初学者的问题会自然消解。

–

Knight, F. H. (1921). Risk, Uncertainty, and Profit. Boston, MA: Hart, Schaffner & Marx.

Tagged Bayesian

2007 11 11 lixiaoxu Chinese Comments (1) Permalink

Understanding the nominal IV

n=10000;r=0.6;r_e=(1-r*r)^.5; X=rnorm(n); Y=X*r+r_e*rnorm(n); Z=as.integer(X>0); Y=Y*Z+(-Y)*(1-Z); Z=as.factor(X>0); ## norminal IV ##Red (totally covered by Green) : Y~X ##Green: Y~X+Z Blue: Y~X+Z+X:Z ##Brown: X~Y+Z Yellow: X~Y+Z+Y:Z plot(X,Y); points(X,predict(lm(Y~X)),col="red"); points(X,predict(lm(Y~X+Z)),col="green"); points(X,predict(lm(Y~X+Z+X:Z)),col="blue"); points(predict(lm(X~Y+Z)),Y,col="brown"); points(predict(lm(X~Y+Z+Y:Z)),Y,col="yellow");

Tagged R, regression

2007 11 09 lixiaoxu English Comments (3) Permalink

答:有同学认为不应该浪费时间教三遍p值和置信区间

如果确实大部分同学认真跟着我学三遍后还不能明白区间估计的假设检验，我承认是我教学上的失败。然而我不介意讲第四遍第五遍（实际上，在结构方程部分， $\Sigma(\theta)$ 、方程结构和S的关系我至少重复了五遍。但是五遍都能听懂，一定胜过三遍还没听懂？）假如有同学有兴趣，欢迎贡献一个问卷调查有多少人终于弄懂区间估计和假设检验，还没有弄懂的同学中有多少同学仍然有足够的兴趣企图花时间去弄懂。做在线问卷只需要动机，不需要写代码的能力。我很希望有更多同学去实践在线问卷这项重要的技能。

到底教什么是重要的，我的判断没有改变。我仍确信选讲p值、区间估计符合我对大家学术倾向的最初预期。倾向学术还是倾向职场，对每位同学无所谓对错。要错就是我最初的预期错。现在确实有同学认为，不搞懂p值、power、区间估计照样可以安心作学术，照样可以面对海量的报告p值的文献，照样可以在自己的学术作品中每篇都报告p值。对这一类同学，我以为这是把学术当作普通谋生行业。我要编量表宁可划这类同学为职场倾向。但如果有同学对p、CI这类貌似非应用的学术问题感兴趣（当然有），我认为太有必要在北大的研究生课程里占用足够的正课时间。这是我的公开立场。

同样，我也相信大家对于什么东西是重要的自有度量。但这并不意味着我的课程需要符合多数同学目前的偏好。我相信许多同学考入北大，是为了有一个机会让北大改变自己的偏好（或品味？），而不是相反。

至于上课提点考试题目的噱头，建议对此有意见的同学把它正确地理解为我的一种调侃方式，而不是我对考题的允诺或背书。同时建议对此调侃方式不能接受的同学，错误地把它理解为我对试题的某种程度的提示，我尽量弄假成真促成喜剧。

最后，我很清楚以上这些颇为偏激的意见显然不适合作为一个comment跟在任何一个同学的学习笔记之后。因为每一位在学习笔记中花费时间陈述课程意见的同学，他们是在为课程作自己的一份义工。义工身后，多的是搭便车的沉默群众。甚至我的comment本身,更多时是基于我对原贴的片面误读而不是全面的解读，因为原贴全篇超过60%的篇幅在正面肯定我的课程教学。显然，我的这篇答复意见已经完全不针对原贴和发表原贴的那位同学，所以我决定把这个回复贴在自己的教学笔记,并欢迎所有同学匿名或者不匿名评论。

Tagged confidence interval, p value

2007 11 02 lixiaoxu Chinese Comments (4) Permalink

R: Simulating multiple normal distribution with any given corr matrix

For example , we have a corr matrix for five standardized factors $\left[\begin{array}{c}1.00\;0.42\;0.41\;0.55\;0.42\\ 0.42\;1.00\;0.48\;0.47\;0.46\\ 0.41\;0.48\;1.00\;0.48\;0.44\\ 0.55\;0.47\;0.48\;1.00\;0.50\\ 0.42\;0.46\;0.44\;0.50\;1.00\end{array}\right]$ (Hau, Chinese Textbook, pp. 49-50).

m_corr=c(1, .42, .41, .55, .42); m_corr=cbind(m_corr, c(.42, 1,.48, .47, .46)); m_corr=cbind(m_corr, c(.41, .48,1, .48, .44)); m_corr=cbind(m_corr, c(.55, .47,.48, 1, .50)); m_corr=cbind(m_corr, c(.42, .46,.44, .50, 1)); m_corr; ## show the original corr matrix dp=svd(m_corr); plot(dp$d,type="o"); ## show the scree plot if there is a PC analysis diag(dp$d); ## show the diagonal matrix of eigen values dp$u %*% diag(dp$d) %*% t(dp$u); ## show corr matrix P=dp$u %*% diag(sqrt(dp$d)); P %*% t(P); ## show corr matrix again n=300; ## typical sample size, you can also try n=30000 f_five=matrix(rnorm(n*5),nrow=n); ## independant multiple normal distribution f_five=f_five %*% t(P); ## correlated cor(f_five); ## test the sample corr matrix

Tagged Normal Distribution, R, Simulation

2007 10 27 lixiaoxu English Comments (1) Permalink

The tail(s) of p value

For any given $H_0$ vs $H_1$ , the p value of any given point x is $\underset{\theta\in H_{0}}{sup}P\left(\left\{ z|L\left(z\right)\ge L\left(x\right)\right\} |\theta\right)$ , Where $L\left(x\right)\equiv\frac{\underset{\theta\in H_{1}}{sup}\left[f\left(x|\theta\right)\right]}{\underset{\theta\in H_{0}}{sup}\left[f\left(x|\theta\right)\right]}$

– See R. Weber’s Statistics Note (Chap 6.2 & 7.1)

I made some wrong comment on the pdf Null Ritual (Gigerenzer, Krauss, & Vitouch, 2004) Where three types of significance level (rather than p value) were discussed. I had written the comment to note that the chapter had ignored the role of $H_1$ in definition of p value. In almost every textbook, the two-tail p vs single-tail p are differentiated. Usually, the two-tail p is defined by $H_1$ like $\mu\neq0$ .

Here I demonstrate a three-tail p value case on R platform.

z=(-1000:1000)*0.02; f=0.5 * dchisq(abs(z),df=5); h=dchisq(10,df=5)*.5; plot(z,f,type="h",col=c("black","grey")[1+(f>h)]); lines(c(-20,20),c(h,h)); ## is * binomial(-1 vs 1) ##

Do you agree the region nearby zero under the “V” curve (which is below the horizontal line) should be the 3rd tail? I think so, if only $H_1$ includes all other possible distributions in the same shape.

You’ll also agree there will be two asymmetrical tails if $H_1$ includes just two asymmetrical curves, for example, $\mu=-2$ and $\mu=1$ ( $\sigma^2\equiv1$ ) while $H_0$ is the standardized normal distribution.