统计软件应用:使用MLwiN作简单两层分析

所用到的数据集popular来源于UCLA统计教学公开网站上HLM教学部分,链接地址为:http://www.ats.ucla.edu/stat/hlm/examples/ma_hox/chapter2.htm

使用技术来源于李晓煦老师心理统计课程教学ppt,模型借鉴了上面网站上的示例(因为使用的是同一组数据)

打开MLwiN

File-> New worksheet

Data Manipulation->view or edit data

Paste

1

Name 修改

2

Model->Equations

3

Graphs

5

Apply

6

接下来用Two-level: The variable sex is included as a random effect and teacher experience (texp) as fixed effect

7

Estimate

8

使用hsb12数据集做简单多层分析

数据为李晓煦老师《心理统计(二)》所提供。

hsb12

取用其中性别(female)和社会地位(ses)对不同学校的学生数学能力作多层分析,学校为第一层变量。分析中拟合了4个不同的嵌套模型,观察不同的拟合效果。

因为要用到nlme软件包,所以在不能直接在Rweb页面运行。现提供代码以及运行结果。

> Dataset<-read.csv(choose.files(),header=T)

> names(Dataset)
> attach(Dataset)

> lme.1<-lme(mathach~1+minority+ses,random=~1+minority+ses|school,data=Dataset)   ##第一个模型不含交叉

> lme.2<-lme(mathach~1+minority+ses+I(ses^2),random=~1+minority+ses|school,data=Dataset)  ##不含交叉但是社会地位包含平方项

> lme.3<-lme(mathach~1+minority+ses+I(ses^2)+I(ses*minority),random=~1+minority+ses|school,data=Dataset)  ##含交叉,社会地位有平方项

> lme.4<-lme(mathach~1+minority+ses+I(ses*minority),random=~1+minority+ses|school,data=Dataset)  ##含交叉,不含平方项

##summary结果如下

> summary(lme.1)
Linear mixed-effects model fit by REML
Data: Dataset
AIC BIC logLik
46449.98 46518.78 -23214.99

Random effects:
Formula: ~1 + minority + ses | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.8605804 (Intr) minrty
minority 1.2425164 0.317
ses 0.5182275 -0.241 -0.829
Residual 5.9828571

Fixed effects: mathach ~ 1 + minority + ses
Value Std.Error DF t-value p-value
(Intercept) 13.491029 0.1746076 7023 77.26487 0
minority -3.076250 0.2422245 7023 -12.70000 0
ses 2.105722 0.1141631 7023 18.44486 0
Correlation:
(Intr) minrty
minority -0.161
ses -0.137 0.021

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.17581712 -0.71576912 0.03386765 0.75304866 3.03705436

Number of Observations: 7185
Number of Groups: 160
> summary(lme.2)
Linear mixed-effects model fit by REML
Data: Dataset
AIC BIC logLik
46454.36 46530.03 -23216.18

Random effects:
Formula: ~1 + minority + ses | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.859616 (Intr) minrty
minority 1.252320 0.317
ses 0.511626 -0.294 -0.863
Residual 5.983493

Fixed effects: mathach ~ 1 + minority + ses + I(ses^2)
Value Std.Error DF t-value p-value
(Intercept) 13.460861 0.18248116 7022 73.76575 0.0000
minority -3.089690 0.24287462 7022 -12.72134 0.0000
ses 2.113992 0.11469128 7022 18.43203 0.0000
I(ses^2) 0.066774 0.09891346 7022 0.67507 0.4997
Correlation:
(Intr) minrty ses
minority -0.140
ses -0.179 0.011
I(ses^2) -0.292 -0.042 0.113

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.18434936 -0.71665805 0.03245726 0.75343379 3.02687521

Number of Observations: 7185
Number of Groups: 160

> summary(lme.3)
Linear mixed-effects model fit by REML
Data: Dataset
AIC BIC logLik
46432.6 46515.15 -23204.3

Random effects:
Formula: ~1 + minority + ses | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.7886997 (Intr) minrty
minority 1.5619471 0.422
ses 0.3249365 -0.482 -0.982
Residual 5.9755974

Fixed effects: mathach ~ 1 + minority + ses + I(ses^2) + I(ses * minority)
Value Std.Error DF t-value p-value
(Intercept) 13.494406 0.1776984 7021 75.93995 0.0000
minority -3.195230 0.2590601 7021 -12.33393 0.0000
ses 2.455783 0.1284693 7021 19.11572 0.0000
I(ses^2) -0.071714 0.1014760 7021 -0.70671 0.4798
I(ses * minority) -1.219800 0.2325023 7021 -5.24640 0.0000
Correlation:
(Intr) minrty ses I(s^2)
minority -0.074
ses -0.149 -0.018
I(ses^2) -0.302 -0.019 -0.049
I(ses * minority) -0.030 0.073 -0.518 0.278

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.21941727 -0.71504662 0.02983522 0.74954540 2.95214218

Number of Observations: 7185
Number of Groups: 160

> summary(lme.4)
Linear mixed-effects model fit by REML
Data: Dataset
AIC BIC logLik
46428.32 46503.99 -23203.16

Random effects:
Formula: ~1 + minority + ses | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.789523 (Intr) minrty
minority 1.559133 0.423
ses 0.349173 -0.501 -0.973
Residual 5.975039

Fixed effects: mathach ~ 1 + minority + ses + I(ses * minority)
Value Std.Error DF t-value p-value
(Intercept) 13.459842 0.1694865 7022 79.41543 0
minority -3.199439 0.2587943 7022 -12.36287 0
ses 2.451657 0.1287472 7022 19.04241 0
I(ses * minority) -1.174167 0.2235735 7022 -5.25181 0
Correlation:
(Intr) minrty ses
minority -0.083
ses -0.181 -0.023
I(ses * minority) 0.060 0.078 -0.525

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.22775400 -0.71472472 0.02825461 0.74934477 2.94494463

Number of Observations: 7185
Number of Groups: 160

##也可以直接观察AIC函数结果,数越小则拟合模型越好。

> AIC(lme.1,lme.4)
df AIC
lme.1 10 46449.98
lme.4 11 46428.32

> AIC(lme.2,lme.4)
df AIC
lme.2 11 46454.36
lme.4 11 46428.32

> AIC(lme.3,lme.4)
df AIC
lme.3 12 46432.60
lme.4 11 46428.32

##可见,最好的模型是

lme.4<-lme(mathach~1+minority+ses+I(ses*minority),random=~1+minority+ses|school,data=Dataset)

《统计软件应用》:用reshape包操纵数据、作图

##相应代码实例在 R in action (Robert I. Kabacoff)书中reshape包章节可见。



《统计软件运用》:绘制中国主要城市地图

In this post, 624 main cited are pointed on the map, which shall indicates the density of cities in China.
Cities from different provinces are pointed with different colors.
And the names of the Provinces are also marked.
The data has been derived from the Internet.

lat_long



运行结果如下:

map

《统计软件应用》作业:编写一个plotEx函数(修改)

作业要求:
写一个plotEx(x,y)函数
用分组密度图表现x,y一个连续一个离散的情形
用气泡图表现x,y两变量都离散的情况
对于x,y都连续的情况,区分
在样本量超过500的时候,用灰度表示区域的点集密度
在样本量未达500的时候,添加参考椭圆、回归线与平滑线(loess)
注:需要car,sm程序包,Rweb上无法显示,需要在Rstudio中运行

plotEx<-function(x,y) {
a<-is.numeric(x)==T
b<-is.numeric(y)==T
n500){
smoothScatter(x,y,main="Scatterplot Colored by Smoothed Densities")
}
else {
library(car)
scatterplot(x,y,ellipse=T,boxplot=F)
abline(lm(y~x))
lines(loess(x,y))
}
}
else if (a & !b){
library(sm)
sm.density.compare(x,y)
}else if (!a & b){
library(sm)
sm.density.compare(y,x)
}
else{
symbols(x,y)
}
}

接下来做几个例子
ChickWeight

ChickWeight包含578个数据.

> data1 attach(data1)
> names(data1)
[1] "weight" "Time" "Chick" "Diet"
> plotEx(Diet, weight)
weight~Diet

> plotEx(Time,weight)
weigh~Time

trees

trees包含31个数据.

> data2 attach(data2)
> names(data2)
[1] "Girth" "Height" "Volume"
> plotEx(Height,Volume)
V~H

> plotEx(Girth,Volume)
V~G

《统计软件应用》:编写一个扩充版summaryEx.data.frame函数

作业要求:
只需要处理factor 与 numeric两种情况
能给出缺失值的个数、非缺失值的个数
能通过unique(…)函数判断一个numeric变量是否在逻辑上是factor变量(取值个数未超过短时记忆上限5),附带给出作为factor变量的概览
对任意一对factor变量报告表格式概览

function 名称:myfun
输入格式:myfun(dataframe$vname)
例子:ChickWeight (from the R Datasets)



《统计应用软件》:cat()函数运用



《统计软件应用》课程作业:试举5例尝试过的R代码中容易出错的情形

1\ 数据分析之前的处理:
#把缺失、无用数据标识为NA:
(1)首先给缺失、无用数据赋值(可用Excel) 如0、99
(2)>data$a[data$a==0]<-NA

#剔除无用数据(注意使用前先将要剔除的范围标识为“NA”):
〉newdata<-na.omit(dataset)

2\ subset的使用(data放在开头,而不像lm()中data放在末尾;以及“|”的使用):
>subset(data,INDEX==1|INDEX==2) ##选择INDEX为1和2的行数据

3\ 分析之前注意数据类型转换(以将numeric指代的factor为例):
dataset$ID<-as.factor(dataset$ID)

4\ 选择观测时,在数据集后用“[]”时,注意不要忘记“,”来分隔选择的与列。比如,在以gender变量来筛选时,要取所有变量进入新数据集,“,”后为空格。
>newdata<-leadership[which(gender=='M'),]

5\ 统一作图时,注意标明x轴范围
>plot(x,y,xlim=c(0,50),ylim=c(0,60))

6\ 图片自动保存
#R的命令行中路径名的使用,注意“/”应为“//”
#jpeg 也可以为 pdf等格式
#注意文件路径中paste()的使用
#不要忘记dev.off()关闭路径
jpeg(paste("C:\\Users\\Administrator\\Desktop\\Graph3\\",cid,".",i,".jpeg",sep=""))
plot(........)
dev.off()

Power of var.test

尝试用图示表示F检验的Power. 假设B组的总体标准差sigma为A组总体标准差的1.5倍. 现从A、B组分别抽样na、nb. 希望展现的是:
1、两组总体的平均值mu1、mu2对结果没有影响;
2、改变na、nb的值,观察检验的效力(Power)变化.

观察独立样本F检验B组sigma为A组1.5倍时的Power
首先设置na=5,nb=10



na=10; nb=10
发现效力增大到30%左右



na=5,nb=100
可以发现效力很低




na=50,nb=50
效力较高
读者可再手动改变na、nb的值,观察效力在90%以上有怎样的规律.
比如na=50,nb=100;na=70,b=nb=70


Han & Chi (2012, p. 1635)

这是李晓煦老师的ppt上的给学生的tasks。数据来自《自发社会比较中的威胁效应及自我平衡策略》(韩晓燕, 迟毓凯)中分析报告4.4.2。
数据内容为:实验任务是否相关(“高相关”和“低相关”)、被试是否经过自我肯定( “有自我肯定”和“无自我肯定”)两个变量对被试的自我评价水平的交叉作用。原始数据只提供了被试总人数,没有提供没有报告各种交叉条件的样本量。根据李晓煦老师的ppt,假设各种交叉条件的样本量为44(高相关、自我肯定)/32(高相关,无自我肯定)/34(低相关,自我肯定)/31(低相关,无自我肯定)

以下摘自原始文献:
“实验任务“高相关”条件下, “无自我肯定”的被试的自我评价水平显著低于“有自我肯定”条件下被试的自我评价(Maffirm=15.94,SD=1.86; Mnoaffirm=13.93, SD=2.75), t=3.80, df=74,p<0.001。在实验任务“低相关”条件下, 虽然“有自我肯定”条件下被试的自我评价高于“无自我肯定”条件下的被试(Maffirm=16.13, SD=2.32; Mnoaffirm=15.78, SD=2.58)”

参考文献:韩晓燕, 迟毓凯,《自发社会比较中的威胁效应及自我平衡策略》,心理学报2012, Vol. 44, No.12, 1628.1640