user — Oct 13, 2014, 10:45 AM
setwd("~/Documents/Universiteit/Psychometrics")
library(foreign)
library(psy)
load("SI.RData")
data <- read.spss("TOTAAL.SAV",to.data.frame=T)
Warning: duplicated levels in factors are deprecated
Warning: duplicated levels in factors are deprecated
Warning: duplicated levels in factors are deprecated
Warning: duplicated levels in factors are deprecated
AsNumeric <- function(var.names,dat=data){
for (j in 1:ncol(dat[var.names])){
dat[var.names][,j] <- as.numeric(dat[var.names][,j])
}
return(dat)
}
#############################
#Disinhibition & Impulsivity#
#############################
#Impulsivity
data <- AsNumeric(Impulsivity87)
data <- AsNumeric(Impulsivity91)
data$S255 <- 8-data$S255
data$S264 <- 8-data$S264
data$S274 <- 8-data$S274
data$S292 <- 8-data$S292
data$SS204 <- 8-data$SS204
data$SS211 <- 8-data$SS211
data$SS219 <- 8-data$SS219
data$SS235 <- 8-data$SS235
data <- AsNumeric(Disinhibition87)
data <- AsNumeric(Disinhibition91)
impulsivity87 <- AsNumeric(Impulsivity87, data)[Impulsivity87]
impulsivity91 <- AsNumeric(Impulsivity91, data)[Impulsivity91]
disinhibition87 <- AsNumeric(Disinhibition87, data)[Disinhibition87]
disinhibition91 <- AsNumeric(Disinhibition91, data)[Disinhibition91]
# Paper 3
require("OpenMx")
Loading required package: OpenMx
require(semPlot)
Loading required package: semPlot
obs.var = c('Dis0', 'Dis1', 'Dis2', 'Dis3', 'Dis4', 'Dis5', 'Imp0', 'Imp1', 'Imp2', 'Imp3', 'Imp4', "Imp5")
#dimnames(cormat) <- list(obs.var, obs.var)
dataMX <- cbind(disinhibition87,impulsivity87)
names(dataMX) <- obs.var
#Create an MxModel object
modelDI <- mxModel("Model",
type="RAM",
#mxData(observed=cormat,type="cor", numObs=500),
mxData(observed=dataMX,type="raw"),
manifestVars=obs.var,
latentVars=c("Disinhibition", "Impulsivity"),
# residual variances
mxPath(from=obs.var,
arrows=2,
free=TRUE,
values=c(1,1,1,1,1,1,1,1,1,1,1,1),
labels=c("e0","e1","e2","e3","e4","e5","e6","e7","e8","e9","e10","e11")
),
# latent variance
mxPath(from="Disinhibition",
arrows=2,
free=FALSE,
values=1,
labels ="varDis"
),
# latent variance
mxPath(from="Impulsivity",
arrows=2,
free=FALSE,
values=1,
labels ="varImp"
),
# latent covariance
mxPath(from="Disinhibition",
to="Impulsivity",
arrows=2,
free=TRUE,
values=1,
labels ="covDI"
),
# factor loadings
mxPath(from="Disinhibition",
to=obs.var[1:6],
arrows=1,
free=c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE),
values=c(1,1,1,1,1,1),
labels =c("l0","l1","l2","l3","l4","l5")
),
# factor loadings
mxPath(from="Impulsivity",
to=obs.var[7:12],
arrows=1,
free=c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE),
values=c(1,1,1,1,1,1),
labels =c("l6","l7","l8","l9","l10","l11")
),
#means and intercepts
mxPath(from="one",
to=c("Disinhibition","Impulsivity"),
arrows=1,
free=TRUE,
values=c(1,1),
labels=c("meanD","meanI")
)
) # close model
fit <- mxRun(modelDI)
Running Model
semPaths(fit, whatLabels="stand")
summary(fit)
data:
$Model.data
Dis0 Dis1 Dis2 Dis3 Dis4
Min. :1.00 Min. :1.0 Min. :1.00 Min. :1.00 Min. :1.00
1st Qu.:2.00 1st Qu.:1.0 1st Qu.:1.00 1st Qu.:5.00 1st Qu.:5.00
Median :4.00 Median :3.0 Median :2.00 Median :6.00 Median :6.00
Mean :3.66 Mean :3.3 Mean :2.71 Mean :5.34 Mean :5.52
3rd Qu.:5.00 3rd Qu.:5.0 3rd Qu.:4.00 3rd Qu.:7.00 3rd Qu.:7.00
Max. :7.00 Max. :7.0 Max. :7.00 Max. :7.00 Max. :7.00
NA's :22 NA's :28 NA's :29 NA's :23 NA's :29
Dis5 Imp0 Imp1 Imp2
Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.00
1st Qu.:5.00 1st Qu.:3.00 1st Qu.:2.00 1st Qu.:3.00
Median :6.00 Median :4.00 Median :4.00 Median :4.00
Mean :5.96 Mean :4.27 Mean :3.72 Mean :4.23
3rd Qu.:7.00 3rd Qu.:6.00 3rd Qu.:5.00 3rd Qu.:5.00
Max. :7.00 Max. :7.00 Max. :7.00 Max. :7.00
NA's :23 NA's :39 NA's :27 NA's :26
Imp3 Imp4 Imp5
Min. :1.00 Min. :1.00 Min. :1.00
1st Qu.:2.00 1st Qu.:1.00 1st Qu.:2.00
Median :4.00 Median :3.00 Median :3.00
Mean :3.59 Mean :3.22 Mean :3.38
3rd Qu.:5.00 3rd Qu.:4.00 3rd Qu.:5.00
Max. :7.00 Max. :7.00 Max. :7.00
NA's :25 NA's :25 NA's :19
free parameters:
name matrix row col Estimate Std.Error
1 l0 A Dis0 Disinhibition 0.7525 0.01921
2 l1 A Dis1 Disinhibition 0.6712 0.01761
3 l2 A Dis2 Disinhibition 0.5494 0.01500
4 l3 A Dis3 Disinhibition 1.0808 0.02518
5 l4 A Dis4 Disinhibition 1.1148 0.02582
6 l5 A Dis5 Disinhibition 1.1853 0.02691
7 l6 A Imp0 Impulsivity 0.8783 0.02566
8 l7 A Imp1 Impulsivity 0.7658 0.02251
9 l8 A Imp2 Impulsivity 0.8641 0.02485
10 l9 A Imp3 Impulsivity 0.7421 0.02195
11 l10 A Imp4 Impulsivity 0.6643 0.02030
12 l11 A Imp5 Impulsivity 0.6972 0.02091
13 e0 S Dis0 Dis0 3.2945 0.11589
14 e1 S Dis1 Dis1 3.4307 0.11971
15 e2 S Dis2 Dis2 3.0692 0.10628
16 e3 S Dis3 Dis3 1.3508 0.05838
17 e4 S Dis4 Dis4 1.1005 0.05182
18 e5 S Dis5 Dis5 1.0377 0.05314
19 e6 S Imp0 Imp0 2.4329 0.09502
20 e7 S Imp1 Imp1 2.0936 0.07985
21 e8 S Imp2 Imp2 2.0261 0.08068
22 e9 S Imp3 Imp3 2.0769 0.07930
23 e10 S Imp4 Imp4 2.9152 0.10484
24 e11 S Imp5 Imp5 2.3803 0.08800
25 covDI S Disinhibition Impulsivity 0.1021 0.03432
26 meanD M 1 Disinhibition 4.9711 0.11651
27 meanI M 1 Impulsivity 4.8593 0.14013
Std.Estimate Std.SE lbound ubound
1 0.3830 0.009775
2 0.3407 0.008938
3 0.2993 0.008169
4 0.6810 0.015866
5 0.7283 0.016867
6 0.7584 0.017221
7 0.4907 0.014336
8 0.4678 0.013750
9 0.5189 0.014923
10 0.4578 0.013542
11 0.3626 0.011083
12 0.4118 0.012351
13 0.8533 0.030019
14 0.8839 0.030842
15 0.9104 0.031527
16 0.5363 0.023177
17 0.4696 0.022115
18 0.4248 0.021757
19 0.7593 0.029654
20 0.7812 0.029795
21 0.7307 0.029099
22 0.7904 0.030180
23 0.8685 0.031236
24 0.8304 0.030703
25 0.1021 0.034320
26 NA NA
27 NA NA
observed statistics: 20985
estimated parameters: 27
degrees of freedom: 20958
-2 log likelihood: 79647
saturated -2 log likelihood: NA
number of observations: 1775
chi-square: NA
p: NA
Information Criteria:
df Penalty Parameters Penalty Sample-Size Adjusted
AIC: 37731 79701 NA
BIC: -77151 79849 79763
CFI: NA
TLI: NA
RMSEA: NA
timestamp: 2014-10-13 10:45:26
frontend time: 0.2469 secs
backend time: 2.293 secs
independent submodels time: 6.986e-05 secs
wall clock time: 2.54 secs
cpu time: 2.54 secs
openmx version number: 1.4-3532
# New paper
library(ltm)
Loading required package: MASS
Loading required package: msm
Loading required package: polycor
Loading required package: mvtnorm
Loading required package: sfsmisc
imp87 <- (impulsivity87>3)*1.0
fit2PL <- ltm(imp87 ~ z1)
plot.ltm(fit2PL, zrange=c(-5,5))
imp87 <- (impulsivity87>4)*1.0
fit2PL <- ltm(imp87 ~ z1)
plot.ltm(fit2PL, zrange=c(-5,5))
# 4 is center
# Now split every question into two
imp87.double <- cbind( impulsivity87>3, impulsivity87>5 ) * 1.0
fit2PL <- ltm(imp87.double ~ z1)
plot.ltm(fit2PL, zrange=c(-5,5))
# Rasch model with unconstrained discrimination parameter
fitRasch <- rasch(imp87.double)
# Two-Parameter Logistic model
fit2PL <- ltm(imp87.double ~ z1)
# Three-Parameter model with discrimination parameter equal to 1
#fit3PL.rc <- tpm(imp87.double, type = "rasch", constraint = cbind(1:5, 3, 1))
# Three-Parameter model with equal discrimination parameter across items
fit3PL.r <- tpm(imp87.double, type = "rasch")
# Three-Parameter model -- no constraint
fit3PL <- tpm(imp87.double)
Warning: Hessian matrix at convergence is not positive definite; unstable solution.
# Anova ltm check
anova.ltm(fit2PL, fit3PL)
Likelihood Ratio Table
AIC BIC log.Lik LRT df p.value
fit2PL 22282 22414 -11117
fit3PL 22304 22501 -11116 2.17 12 0.999
tcc<-function(fit){# begin Test Characteristic Function
n<-dim(fit$X)[2]; p<-plot(fit, items=1, type="ICC", plot=FALSE)
r<-length(p[,2]); X<-p[,1]; iccs<-matrix(NA,ncol=n,nrow=r)
for (i in 1:n) iccs[,i]<-plot(fit,items=i,type="ICC", plot=FALSE)[,2]
Y<-apply(iccs, MARGIN=1, sum);
cbind(X,Y)
} # end Test Characteristic Function
plot(tcc(fitRasch))
plot(tcc(fit2PL))
plot(tcc(fit3PL))
plot.ltm(fit2PL, type="IIC", zrange=c(-6,6))
tiic<-function(fit){
n<-dim(fit$X)[2]; p<-plot(fit, items=1, type="IIC", plot=FALSE)
r<-length(p[,2]); X<-p[,1]; iccs<-matrix(NA,ncol=n,nrow=r)
for (i in 1:n) iccs[,i]<-plot(fit,items=i,type="IIC", plot=FALSE)[,2]
Y<-apply(iccs, MARGIN=1, sum);
cbind(X,Y)
}
plot(tiic(fit2PL), type="l")
require("OpenMx")
fit <- fit2PL
n<-dim(fit$X)[2]; p<-plot(fit, items=1, type="ICC", plot=FALSE)
r<-length(p[,2]); X<-p[,1]; iccs<-matrix(NA,ncol=n,nrow=r)
for (i in 1:n) iccs[,i]<-plot(fit,items=i,type="ICC", plot=FALSE)[,2]
Y<-apply(iccs, MARGIN=1, sum);
#Create an MxModel object
oneFactorModel <- mxModel("Common Factor Model Path Specification",
type="RAM",
mxData(observed=LSAT,type="raw"),
manifestVars=c("Item 1", "Item 2", "Item 3", "Item 4", "Item 5"),
latentVars="Factor",
# residual variances
mxPath(from=c("Item 1", "Item 2", "Item 3", "Item 4", "Item 5"),
arrows=2,
free=TRUE,
values=c(1,1,1,1,1),
labels=c("e1","e2","e3","e4","e5")
),
# latent variance
mxPath(from="Factor",
arrows=2,
free=FALSE,
values=1,
labels ="varFactor"
),
# factor loadings
mxPath(from="Factor",
to=c("Item 1", "Item 2", "Item 3", "Item 4", "Item 5"),
arrows=1,
free=c(TRUE,TRUE,TRUE,TRUE,TRUE),
values=c(1,1,1,1,1),
labels =c("l1","l2","l3","l4","l5")
),
# means
mxPath(from="one",
to=c("Item 1", "Item 2", "Item 3", "Item 4", "Item 5","Factor"),
arrows=1,
free=c(TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),
values=c(1,1,1,1,1,0),
labels =c("meanx1","meanx2","meanx3","meanx4","meanx5",NA)
)
) # close model
oneFactorFit <- mxRun(oneFactorModel)
Running Common Factor Model Path Specification
estimates<-oneFactorFit@output$estimate
lambda<-matrix(estimates[1:5],nrow=1)
mu<-matrix(estimates[11:15],nrow=1)
X<-matrix(X,ncol=1)
ones<-matrix(1,ncol=1,nrow=dim(X)[1])
# transform
expected.responses<-ones%*%mu+ X%*%lambda
Y.FA<-apply(expected.responses,1,sum)
plot(X,Y)