Dear list members,
Could any expert on factor analysis be so kind to explain how to calculate AIC
on the output of factanal. Do I calculate AIC wrong or is
factanal$criteria["objective"] not a negative log-likelihood?
Best regards
Jens Oehlschl?gel
The AIC calculated using summary.factanal below don't appear correct to me:
n items factors total.df rest.df model.df LL AIC AICc
BIC
1 1000 20 1 210 170 40 -5.192975386 90.38595 93.80618
286.6962
2 1000 20 2 210 151 59 -3.474172303 124.94834 132.48026
414.5059
3 1000 20 3 210 133 77 -1.745821627 157.49164 170.51984
535.3888
4 1000 20 4 210 116 94 -0.120505369 188.24101 207.97582
649.5700
5 1000 20 5 210 100 110 -0.099209921 220.19842 247.66749
760.0515
6 1000 20 6 210 85 125 -0.072272695 250.14455 286.18574
863.6140
7 1000 20 7 210 71 139 -0.054668588 278.10934 323.36515
960.2873
8 1000 20 8 210 58 152 -0.041708051 304.08342 358.99723
1050.0622
9 1000 20 9 210 46 164 -0.028686298 328.05737 392.87174
1132.9292
10 1000 20 10 210 35 175 -0.015742783 350.03149 424.78877
1208.8887
11 1000 20 11 210 25 185 -0.007007901 370.01402 454.55947
1277.9487
12 1000 20 12 210 16 194 -0.005090725 388.01018 481.99776
1340.1147
summary.factanal <- function(object, ...){
if (inherits(object, "try-error")){
c(n=NA, items=NA, factors=NA, total.df=NA, rest.df=NA, model.df=NA, LL=NA,
AIC=NA, AICc=NA, BIC=NA)
}else{
n <- object$n.obs
p <- length(object$uniquenesses)
m <- object$factors
model.df <- (p*m) + (m*(m+1))/2 + p - m^2
total.df <- p*(p+1)/2
rest.df <- total.df - model.df # = object$dof
LL <- -as.vector(object$criteria["objective"])
k <- model.df
aic <- 2*k - 2*LL
aicc <- aic + (2*k*(k+1))/(n-k-1)
bic <- k*log(n) - 2*LL
c(n=n, items=p, factors=m, total.df=total.df, rest.df=rest.df,
model.df=model.df, LL=LL, AIC=aic, AICc=aicc, BIC=bic)
}
}
multifactanal <- function(factors=1:3, ...){
names(factors) <- factors
ret <- lapply(factors, function(factors){
try(factanal(factors=factors, ...))
})
class(ret) <- "multifactanal"
ret
}
summary.multifactanal <- function(object,...){
do.call("rbind", lapply(object, summary.factanal))
}
print.multifactanal <- function(x,...){
ret <- summary.multifactanal(x)
print(ret, ...)
invisible(ret)
}
# simulate a true 4-factor model
n <- 1000
ktrue <- 4
kfac <- 5
true <- matrix(rnorm(n*ktrue), ncol=ktrue)
x <- matrix(rep(true, kfac)+rnorm(n*ktrue*kfac), ncol=ktrue*kfac)
dimnames(x) <- list(NULL, paste(rep(letters[1:ktrue], kfac), rep(1:kfac,
rep(ktrue, kfac)), sep=""))
covmat <- cov.wt(x)
# run factanal for several numbers of factors
mf <- multifactanal(factors=1:12, covmat=covmat)
mf
version
_
platform i386-pc-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 2
minor 5.0
year 2007
month 04
day 23
svn rev 41293
language R
version.string R version 2.5.0 (2007-04-23)
--
"Feel free" - 10 GB Mailbox, 100 FreeSMS/Monat ...
Daniel, Thanks for answering.> your AIC is monotonic increasingThat was the reason I suspected something is wrong, but what? Either the way I calculate the number of estimiated parameters? Or the (scale of the) Likelihood itself?> Have you tried a dimensional reduction technique or to investigate the wieghtings on your paramaters (eigenvalues in spectral decomp - principalcomponents analysis). Maybe some of those factors are not contributing? The data was simulated, so we know what it is. Best regards Jens Oehlschl?gel -- "Feel free" - 10 GB Mailbox, 100 FreeSMS/Monat ...
Dear R Developers:
What is the proper log(likelihood) for 'factanal'? I believe it
should be something like the following:
(-n/2)*(k*(log(2*pi)+1)+log(det(S)))
or without k*(log(2*pi)-1):
(-n/2)*log(det(S)),
where n = the number of (multivariate) observations.
The 'factanal' help file say the factanal component
"criteria:
The results of the optimization: the value of the negative
log-likelihood and information on the iterations used."
However, I'm not able to get this. If it's a log(likelihood),
then replacing a data set m1 by rbind(m1, m1) should double the
log(likelihood). However it has no impact on it. Consider the
following modification of the first example in the 'factanal' help page:
> v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)
> v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
> v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)
> v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
> v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)
> v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
> m1 <- cbind(v1,v2,v3,v4,v5,v6)
> tst <- factanal(m1, factors=3)
> fit1 <- factanal(m1, factors=3)
> fit2 <- factanal(rbind(m1, m1), factors=3)
> fit1$criteria["objective"]
objective
0.4755156
> fit2$criteria["objective"]
objective
0.4755156
From the following example, it seems that the "k*(log(2*pi)-1)
term is omitted:
> x2 <- c(-1,1)
> X2.4 <- as.matrix(expand.grid(x2, x2, x2, x2))
> factanal(X2.4, 1)$criteria
objective counts.function counts.gradient
0 7 7
However, I can't get the 'fit1$criteria' above, even ignoring
the
sample size. Consider the following:
> S3 <- tcrossprod(fit1$loadings)+diag(fit1$uniquenesses)
> log(det(S3))
[1] -6.725685
Shouldn't this be something closer to the 0.4755 for fit1$criteria?
Thanks,
Spencer Graves
JENS: For your purposes, I suggest you try to compute
(-n/2)*log(det(tcrossprod(fit$loadings)+diag(fit$uniquenesses)). If you
do this, you might get more sensible results with your examples.
Hope this helps.
Spencer Graves
Jens Oehlschl?gel wrote:> Dear list members,
>
> Could any expert on factor analysis be so kind to explain how to calculate
AIC on the output of factanal. Do I calculate AIC wrong or is
factanal$criteria["objective"] not a negative log-likelihood?
>
> Best regards
>
>
> Jens Oehlschl?gel
>
>
>
> The AIC calculated using summary.factanal below don't appear correct to
me:
>
> n items factors total.df rest.df model.df LL AIC
AICc BIC
> 1 1000 20 1 210 170 40 -5.192975386 90.38595
93.80618 286.6962
> 2 1000 20 2 210 151 59 -3.474172303 124.94834
132.48026 414.5059
> 3 1000 20 3 210 133 77 -1.745821627 157.49164
170.51984 535.3888
> 4 1000 20 4 210 116 94 -0.120505369 188.24101
207.97582 649.5700
> 5 1000 20 5 210 100 110 -0.099209921 220.19842
247.66749 760.0515
> 6 1000 20 6 210 85 125 -0.072272695 250.14455
286.18574 863.6140
> 7 1000 20 7 210 71 139 -0.054668588 278.10934
323.36515 960.2873
> 8 1000 20 8 210 58 152 -0.041708051 304.08342
358.99723 1050.0622
> 9 1000 20 9 210 46 164 -0.028686298 328.05737
392.87174 1132.9292
> 10 1000 20 10 210 35 175 -0.015742783 350.03149
424.78877 1208.8887
> 11 1000 20 11 210 25 185 -0.007007901 370.01402
454.55947 1277.9487
> 12 1000 20 12 210 16 194 -0.005090725 388.01018
481.99776 1340.1147
>
>
> summary.factanal <- function(object, ...){
> if (inherits(object, "try-error")){
> c(n=NA, items=NA, factors=NA, total.df=NA, rest.df=NA, model.df=NA,
LL=NA, AIC=NA, AICc=NA, BIC=NA)
> }else{
> n <- object$n.obs
> p <- length(object$uniquenesses)
> m <- object$factors
> model.df <- (p*m) + (m*(m+1))/2 + p - m^2
> total.df <- p*(p+1)/2
> rest.df <- total.df - model.df # = object$dof
> LL <- -as.vector(object$criteria["objective"])
> k <- model.df
> aic <- 2*k - 2*LL
> aicc <- aic + (2*k*(k+1))/(n-k-1)
> bic <- k*log(n) - 2*LL
> c(n=n, items=p, factors=m, total.df=total.df, rest.df=rest.df,
model.df=model.df, LL=LL, AIC=aic, AICc=aicc, BIC=bic)
> }
> }
>
> multifactanal <- function(factors=1:3, ...){
> names(factors) <- factors
> ret <- lapply(factors, function(factors){
> try(factanal(factors=factors, ...))
> })
> class(ret) <- "multifactanal"
> ret
> }
>
> summary.multifactanal <- function(object,...){
> do.call("rbind", lapply(object, summary.factanal))
> }
>
> print.multifactanal <- function(x,...){
> ret <- summary.multifactanal(x)
> print(ret, ...)
> invisible(ret)
> }
>
> # simulate a true 4-factor model
> n <- 1000
> ktrue <- 4
> kfac <- 5
> true <- matrix(rnorm(n*ktrue), ncol=ktrue)
> x <- matrix(rep(true, kfac)+rnorm(n*ktrue*kfac), ncol=ktrue*kfac)
> dimnames(x) <- list(NULL, paste(rep(letters[1:ktrue], kfac), rep(1:kfac,
rep(ktrue, kfac)), sep=""))
> covmat <- cov.wt(x)
> # run factanal for several numbers of factors
> mf <- multifactanal(factors=1:12, covmat=covmat)
> mf
>
>
> version
> _
> platform i386-pc-mingw32
> arch i386
> os mingw32
> system i386, mingw32
> status
> major 2
> minor 5.0
> year 2007
> month 04
> day 23
> svn rev 41293
> language R
> version.string R version 2.5.0 (2007-04-23)
>
The log likelihood is
const + n/2* [ log(det(Sigma^-1)) - trace(Sigma^-1*S) ] where Sigma is the
population covariance matrix
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at
stat.math.ethz.ch] On Behalf Of Spencer Graves
Sent: Friday, May 04, 2007 9:20 PM
To: R-devel mailing list
Cc: Jens Oehlschl?gel; r-help at stat.math.ethz.ch
Subject: Re: [R] factanal AIC?
Dear R Developers:
What is the proper log(likelihood) for 'factanal'? I believe it
should be something like the following:
(-n/2)*(k*(log(2*pi)+1)+log(det(S)))
or without k*(log(2*pi)-1):
(-n/2)*log(det(S)),
where n = the number of (multivariate) observations.
The 'factanal' help file say the factanal component
"criteria:
The results of the optimization: the value of the negative log-likelihood and
information on the iterations used."
However, I'm not able to get this. If it's a log(likelihood),
then replacing a data set m1 by rbind(m1, m1) should double the log(likelihood).
However it has no impact on it. Consider the following modification of the
first example in the 'factanal' help page:
> v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)
> v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)
> v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)
> v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)
> v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)
> v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)
> m1 <- cbind(v1,v2,v3,v4,v5,v6)
> tst <- factanal(m1, factors=3)
> fit1 <- factanal(m1, factors=3)
> fit2 <- factanal(rbind(m1, m1), factors=3) >
fit1$criteria["objective"] objective
0.4755156
> fit2$criteria["objective"]
objective
0.4755156
From the following example, it seems that the "k*(log(2*pi)-1) term
is omitted:
> x2 <- c(-1,1)
> X2.4 <- as.matrix(expand.grid(x2, x2, x2, x2)) > factanal(X2.4,
1)$criteria
objective counts.function counts.gradient
0 7 7
However, I can't get the 'fit1$criteria' above, even ignoring
the sample size. Consider the following:
> S3 <- tcrossprod(fit1$loadings)+diag(fit1$uniquenesses)
> log(det(S3))
[1] -6.725685
Shouldn't this be something closer to the 0.4755 for fit1$criteria?
Thanks,
Spencer Graves
JENS: For your purposes, I suggest you try to compute
(-n/2)*log(det(tcrossprod(fit$loadings)+diag(fit$uniquenesses)). If you do
this, you might get more sensible results with your examples.
Hope this helps.
Spencer Graves
Jens Oehlschl?gel wrote:> Dear list members,
>
> Could any expert on factor analysis be so kind to explain how to calculate
AIC on the output of factanal. Do I calculate AIC wrong or is
factanal$criteria["objective"] not a negative log-likelihood?
>
> Best regards
>
>
> Jens Oehlschl?gel
>
>
>
> The AIC calculated using summary.factanal below don't appear correct to
me:
>
> n items factors total.df rest.df model.df LL AIC
AICc BIC
> 1 1000 20 1 210 170 40 -5.192975386 90.38595
93.80618 286.6962
> 2 1000 20 2 210 151 59 -3.474172303 124.94834
132.48026 414.5059
> 3 1000 20 3 210 133 77 -1.745821627 157.49164
170.51984 535.3888
> 4 1000 20 4 210 116 94 -0.120505369 188.24101
207.97582 649.5700
> 5 1000 20 5 210 100 110 -0.099209921 220.19842
247.66749 760.0515
> 6 1000 20 6 210 85 125 -0.072272695 250.14455
286.18574 863.6140
> 7 1000 20 7 210 71 139 -0.054668588 278.10934
323.36515 960.2873
> 8 1000 20 8 210 58 152 -0.041708051 304.08342
358.99723 1050.0622
> 9 1000 20 9 210 46 164 -0.028686298 328.05737
392.87174 1132.9292
> 10 1000 20 10 210 35 175 -0.015742783 350.03149
424.78877 1208.8887
> 11 1000 20 11 210 25 185 -0.007007901 370.01402
454.55947 1277.9487
> 12 1000 20 12 210 16 194 -0.005090725 388.01018
481.99776 1340.1147
>
>
> summary.factanal <- function(object, ...){
> if (inherits(object, "try-error")){
> c(n=NA, items=NA, factors=NA, total.df=NA, rest.df=NA, model.df=NA,
LL=NA, AIC=NA, AICc=NA, BIC=NA)
> }else{
> n <- object$n.obs
> p <- length(object$uniquenesses)
> m <- object$factors
> model.df <- (p*m) + (m*(m+1))/2 + p - m^2
> total.df <- p*(p+1)/2
> rest.df <- total.df - model.df # = object$dof
> LL <- -as.vector(object$criteria["objective"])
> k <- model.df
> aic <- 2*k - 2*LL
> aicc <- aic + (2*k*(k+1))/(n-k-1)
> bic <- k*log(n) - 2*LL
> c(n=n, items=p, factors=m, total.df=total.df, rest.df=rest.df,
model.df=model.df, LL=LL, AIC=aic, AICc=aicc, BIC=bic)
> }
> }
>
> multifactanal <- function(factors=1:3, ...){
> names(factors) <- factors
> ret <- lapply(factors, function(factors){
> try(factanal(factors=factors, ...))
> })
> class(ret) <- "multifactanal"
> ret
> }
>
> summary.multifactanal <- function(object,...){
> do.call("rbind", lapply(object, summary.factanal)) }
>
> print.multifactanal <- function(x,...){
> ret <- summary.multifactanal(x)
> print(ret, ...)
> invisible(ret)
> }
>
> # simulate a true 4-factor model
> n <- 1000
> ktrue <- 4
> kfac <- 5
> true <- matrix(rnorm(n*ktrue), ncol=ktrue) x <- matrix(rep(true,
> kfac)+rnorm(n*ktrue*kfac), ncol=ktrue*kfac)
> dimnames(x) <- list(NULL, paste(rep(letters[1:ktrue], kfac),
> rep(1:kfac, rep(ktrue, kfac)), sep="")) covmat <- cov.wt(x) #
run
> factanal for several numbers of factors mf <-
> multifactanal(factors=1:12, covmat=covmat) mf
>
>
> version
> _
> platform i386-pc-mingw32
> arch i386
> os mingw32
> system i386, mingw32
> status
> major 2
> minor 5.0
> year 2007
> month 04
> day 23
> svn rev 41293
> language R
> version.string R version 2.5.0 (2007-04-23)
>
______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.