https://github.com/cran/Rela
Tip revision: bf9486d1a7621e375b8b8b217bf4483442218699 authored by Michael Chajewski on 25 October 2009, 00:00:00 UTC
version 4.1
version 4.1
Tip revision: bf9486d
itemanal.R
itemanal <- function(object, SE.par=1.96, boots=1000){
if (is.null(row.names(t(object)))) {print("Dataset must be coerced into matrix from prior data frame.")}
else {
if (!is.numeric(object)) {print("Your dataset is not a numeric object.")}
else {
object<-(na.omit(object))
yy<-sort(abs(as.numeric(var(object))))
yy<-yy[1]
yy<-ifelse(yy==0, NA, yy)
if (is.na(yy)) {print("One of your variables is a constant. Constants are disallowed as part of a scale.")}
else {
if ( dim(cbind(object))[2]< 2 ) {print("Your set contains only one variable. At least two are required.")}
else {
if ( dim(cbind(object))[1]< 2 ) {print("Your set contains only one case. You need at least two cases.")}
else {
same <- function(x) {
for (i in 1:ncol(x)) {
for (j in 1:ncol(x)) {
if (i!=j) {
test <- x[,i]-x[,j]
if (sum(abs(test))==0) {print(paste("WARNING: Items ",i," (",colnames(x)[i],") ","and ",j," (",colnames(x)[j],") ","are duplicates."))}
}
}
}
test
}
notneed <- same(object)
object<-(na.omit(object))
options(digits=5)
Type<-mode(object)
Cases<-apply(object,2,length)
Minimum<-apply(object, 2, min)
Maximum<-apply(object, 2, max)
Mean<-apply(object, 2, mean)
Median<-apply(object, 2, median)
sum1<-sum((object-mean(object))^3)
N1<-(length(object))-1
SDF<- function(object) {
sqrt((sum((object-mean(object))^2))/N1)}
SD <- apply(object, 2, SDF)
Sum<-apply(object,2,sum)
SE.mean<-(SD)/(sqrt(length(object)))
Lower<-(Mean)-(SE.par*SD)
Upper<-(Mean)+(SE.par*SD)
Skew<- function(object){
(sum1/N1)*((SDF(object))^3)}
Skewness<-apply(object, 2, Skew)
er.skew <- sqrt(6/(dim(object)[1]))
lower.skew<- SE.par*er.skew
upper.skew<- -SE.par*er.skew
sum2<-sum((object-mean(object))^4)
Kurt<- function(object) {
(sum2/N1)*((SDF(object))^4)}
Kurtosis<-apply(object, 2, Kurt)
er.kurt <- sqrt(24/(dim(object)[1]))
lower.kurt<- SE.par*er.kurt
upper.kurt<- -SE.par*er.kurt
eivar<- apply(object, 2, var)
summary.statement1<-data.frame(cbind("Variable.type"=Type, "No.of.cases"=Cases, Minimum, Maximum, Sum))
summary.statement2<-cbind(Mean, Median, SD, SE.mean, Lower, Upper, "Variance"=eivar)
summary.statement3a<-cbind(Skewness,"SE.skew"=er.skew, "Lower"=lower.skew, "Upper"=upper.skew)
summary.statement3b<-cbind(Kurtosis,"SE.kurtosis"=er.kurt, "Lower"=lower.kurt, "Upper"=upper.kurt)
Covariances<-cov(object,object)
Correlations<-cor(object,object)
# Cronbach's Alpha
Columns<-ncol(object)
Rows<-nrow(object)
Alpha<-(Columns/(Columns-1))*(1-sum(apply(object,2,var))/(var(apply(object,1,sum))))
# Duhachek SE
numbitem = Columns
numbsubj = Rows
itemcov = cov(object)
one=matrix(1,numbitem,1)
jtphij = t(one)
jtphij=jtphij%*%itemcov
jtphij=jtphij%*%one
trmy=sum(diag(itemcov))
trmy=trmy/jtphij
myalpha=1-trmy
nn1=numbitem-1
nn1=numbitem/nn1
myalpha<-nn1*myalpha
trphisq=itemcov%*%itemcov
trphisq=sum(diag(trphisq))
trsqphi=sum(diag(itemcov))
trsqphi=trsqphi^2
ttp=itemcov%*%itemcov
jtphisqj=t(one)
jtphisqj=jtphisqj%*%ttp
jtphisqj=jtphisqj%*%one
omega=trphisq+trsqphi
omega=jtphij*omega
omegab=sum(diag(itemcov))
omegab=omegab*jtphisqj
omega=omega-(2*omegab)
omega=(2/(jtphij^3))*omega
s2=(numbitem^2)/((numbitem-1)^2)
s2=s2*omega
se<-sqrt(s2/numbsubj)
cimin95<-myalpha-(SE.par*se)
cimax95<-myalpha+(SE.par*se)
summary.statement4<-data.frame("No.of.items"=Columns, "Cronbachs.alpha"=Alpha)
summary.statement5<-data.frame("SE.for.alpha"=se, "Lower"=cimin95, "Upper"=cimax95)
# Duhachek SE standardized alpha
itemcov = cor(object)
one=matrix(1,numbitem,1)
jtphij = t(one)
jtphij=jtphij%*%itemcov
jtphij=jtphij%*%one
trmy=sum(diag(itemcov))
trmy=trmy/jtphij
myalpha=1-trmy
nn1=numbitem-1
nn1=numbitem/nn1
myalpha<-nn1*myalpha
trphisq=itemcov%*%itemcov
trphisq=sum(diag(trphisq))
trsqphi=sum(diag(itemcov))
trsqphi=trsqphi^2
ttp=itemcov%*%itemcov
jtphisqj=t(one)
jtphisqj=jtphisqj%*%ttp
jtphisqj=jtphisqj%*%one
omega=trphisq+trsqphi
omega=jtphij*omega
omegab=sum(diag(itemcov))
omegab=omegab*jtphisqj
omega=omega-(2*omegab)
omega=(2/(jtphij^3))*omega
s2=(numbitem^2)/((numbitem-1)^2)
s2=s2*omega
se<-sqrt(s2/numbsubj)
cimin95<-myalpha-(SE.par*se)
cimax95<-myalpha+(SE.par*se)
summary.statement6<-data.frame("No.of.items"=Columns, "Standardized.alpha"=myalpha)
summary.statement7<-data.frame("SE.for.standardized.alpha"=se, "Lower"=cimin95, "Upper"=cimax95)
if.deleted <- function(object) {
item.sum <- apply(object, 2, sum)
scale.sum <- sum(item.sum)
scale.mean <- (scale.sum)/(length(apply(object,1,length)))
scale.mean.var <- data.frame("Variable"=row.names(t(object)), scale.mean)
n.0 <- length(apply(object,2,length))
n.10 <- length(apply(object,1,length))
ave <- function(x) {
(scale.sum-x)/(length(apply(object,1,length)))}
mean.if.deleted <- c()
for (i in 1:length(item.sum)) {
mean.if.deleted<-rbind(mean.if.deleted,ave(item.sum[i]))
}
deleted <- data.frame("Scale.mean"=scale.mean.var[,-1],"Mean.if.deleted"=mean.if.deleted[,1])
}
summary.statement8a <- if.deleted(object)
# Variance If Item Deleted
if.deleted.variance <- function(object) {
item.sum.variance <- apply(object, 1, sum)
item.variance <- var(item.sum.variance)
if (dim(object)[2]<3) {del.variance <- c(var((object)[,1]),var((object)[,2]))
item.variance <- c(item.variance, item.variance)
del.variance <- data.frame("Scale.variance"=item.variance, "Variance.if.deleted"=del.variance)
}
else {
del.variance <- c(1)
for (i in 1:dim(object)[2]) {del.variance<-(c(del.variance,(var(apply(object[,-1*i],1,sum)))))}
del.variance <- del.variance[-1]
del.variance <- data.frame("Scale.variance"=item.variance, "Variance.if.deleted"=del.variance)
}
}
summary.statement8b <- if.deleted.variance(object)
# Corrected total item correlation if item deleted
item.cor<- function(object) {
if (dim(object)[2]<3) {
cor1<- cor(object[,1], object[,2])
cor2<- cor(object[,2], object[,1])
del.cor<- c(cor1, cor2)
del.cor<-data.frame("Corrected.item.total.cor"=del.cor)
}
else
{
del.cor <- c(1)
for (i in 1:dim(object)[2]){ del.cor <- c(del.cor, cor((sapply(object[,i],sum)), (apply(object[,-1*i],1,sum))))}
del.cor <- del.cor[-1]
del.cor<- data.frame("Corrected.item.total.cor"=del.cor)
}
}
summary.statement8 <- data.frame(summary.statement8a, summary.statement8b)
summary.statement9a <- data.frame("Variable"=row.names(t(object)), item.cor(object))
# Square Multiple correlation If Item Deleted
r.if.del<- function(object) {
if (dim(object)[2]<3) {
del.lm<-lm(object[,1]~object[,2])
del.lm.sum1<-summary(del.lm)
del.lm<-lm(object[,2]~object[,1])
del.lm.sum2<-summary(del.lm)
r.if.item.del<-c(del.lm.sum1$r.squared, del.lm.sum2$r.squared)
}
else {
r.if.item.del<- c(1)
for (i in 1:dim(object)[2]) {
r.if.item.del<- c(r.if.item.del, (summary(lm ( object[,i] ~ object[,-1*i])))$r.squared ) }
r.if.item.del <- r.if.item.del[-1]
}
}
summary.statement10<-data.frame("Squared.multiple.cor"=r.if.del(object))
# Alpha if item deleted
del.alpha <- function(object) {
if (dim(object)[2]<3) {make.alpha<-c("N/A","N/A")
}
else {
make.alpha <- c(1)
n.items <- dim(object)[2]
for (i in 1:n.items) {
make.alpha <- c(make.alpha ,(((length(1:n.items)-1) / ((length(1:n.items))-2)) * (1-sum(apply(object[,-1*i],2,var)) / (var(apply(object[,-1*i],1,sum))))))
}
make.alpha <- make.alpha[-1]
}
}
summary.statement9b <- del.alpha(object)
summary.statement9 <- data.frame(summary.statement9a,summary.statement10, "Alpha.if.deleted"=summary.statement9b)
# Bootstrap based alpha standard error and confidence interval
boot.one <-
function(object, rep=boots, ci.par=SE.par, nd=4) {
new <- c(1)
func <- function(object) {
Columns<-ncol(object)
Rows<-nrow(object)
numbitem = Columns
numbsubj = Rows
itemcov = cov(object)
one=matrix(1,numbitem,1)
jtphij = t(one)
jtphij=jtphij%*%itemcov
jtphij=jtphij%*%one
trmy=sum(diag(itemcov))
trmy=trmy/jtphij
myalpha=1-trmy
nn1=numbitem-1
nn1=numbitem/nn1
myalpha<-nn1*myalpha
myalpha }
{ repeat
{ord <- sample(1:length(object[,1]), length(object[,1]), replace=TRUE )
XX <- object[ord,]
XXz <- round(func(XX), digits=nd)
new <- c(new, XXz)
if(length(new)+.1>rep+1){break} }}
new <- new[-1]
orig <- round(func(object), digits=nd)
SE <- round(sd(new), digits=nd)
ci.lower <- round(orig-(ci.par*SE), digits=nd)
ci.upper <- round(orig+(ci.par*SE), digits=nd)
MM <- round(mean(new), digits=nd)
final <- list("new"=new,MM,SE, ci.lower, ci.upper)
final}
b.o.o <- boot.one(object)
b.o.o.sim <- as.numeric(b.o.o$new)
summary.statement10a <- data.frame(cbind(b.o.o[2],b.o.o[3],b.o.o[4],b.o.o[5]))
colnames(summary.statement10a) <- c("Mean", "SE", "Lower", "Upper")
# Bootstrap based standardized alpha standard error and confidence interval
boot.two <-
function(object, rep=boots, ci.par=SE.par, nd=4) {
new <- c(1)
func <- function(object) {
Columns<-ncol(object)
Rows<-nrow(object)
numbitem = Columns
numbsubj = Rows
itemcov = cor(object)
one=matrix(1,numbitem,1)
jtphij = t(one)
jtphij=jtphij%*%itemcov
jtphij=jtphij%*%one
trmy=sum(diag(itemcov))
trmy=trmy/jtphij
myalpha=1-trmy
nn1=numbitem-1
nn1=numbitem/nn1
myalpha<-nn1*myalpha
myalpha }
{ repeat
{ord <- sample(1:length(object[,1]), length(object[,1]), replace=TRUE )
XX <- object[ord,]
XXz <- round(func(XX), digits=nd)
new <- c(new, XXz)
if(length(new)+.1>rep+1){break} }}
new <- new[-1]
orig <- round(func(object), digits=nd)
SE <- round(sd(new), digits=nd)
ci.lower <- round(orig-(ci.par*SE), digits=nd)
ci.upper <- round(orig+(ci.par*SE), digits=nd)
MM <- round(mean(new), digits=nd)
final <- list("new"=new,MM,SE, ci.lower, ci.upper)
final}
b.o.o2 <- boot.two(object)
b.o.o2.sim <- as.numeric(b.o.o2$new)
summary.statement11 <- data.frame(cbind(b.o.o2[2],b.o.o2[3],b.o.o2[4],b.o.o2[5]))
colnames(summary.statement11) <- c("Mean", "SE", "Lower", "Upper")
item.analysis.output<-list(Variables=summary.statement1, Tendency=summary.statement2,
Skewness=summary.statement3a, Kurtosis=summary.statement3b, Covariance=Covariances, Correlation=Correlations,
Alpha=summary.statement4, Conf.Alpha=summary.statement5,Bootstrap.Simmulations=b.o.o.sim,
Alpha.Bootstrap=summary.statement10a, Std.Alpha=summary.statement6, Conf.Std.Alpha=summary.statement7,
Bootstrap.Std.Simmulations=b.o.o2.sim,Std.Alpha.Bootstrap=summary.statement11,
Scale.Stats=summary.statement8, Alpha.Stats= summary.statement9)
ob <- as.character(match.call()[2])
cl <- call("itemanal","object"=ob,"SE.par"=SE.par, "boots"=boots)
item.analysis.output$call <- cl
item.analysis.output$items <- names(item.analysis.output)
class(item.analysis.output) <- c("itemanal", class(item.analysis.output))
item.analysis.output
}}}}}}