args <- commandArgs(trailingOnly=TRUE); input_file <- if(length(args)>=1) args[1] else "bfi-original.csv"
items <- c(paste0("A",1:5),paste0("C",1:5),paste0("E",1:5),paste0("N",1:5),paste0("O",1:5)); rev_items <- c("A1","C4","C5","E1","E2","O2","O5")
runtime_md5 <- function(p) unname(tools::md5sum(p)); source_sha256 <- "68ae71a96c2157b0c49b8d8f8ccaee1292f5e92adce9a3f202fdc7efac883dc9"; analysis_seed <- 20260916L
raw <- read.csv(input_file,check.names=FALSE,na.strings=c("","NA")); stopifnot(all(items%in%names(raw)))
xr <- raw[,items,drop=FALSE]; xr[] <- lapply(xr,as.numeric); z <- unlist(xr); stopifnot(all(z[!is.na(z)]>=1 & z[!is.na(z)]<=6))
n0 <- nrow(xr); miss <- colSums(is.na(xr)); keep <- complete.cases(xr); x <- xr[keep,,drop=FALSE]; n <- nrow(x); nex <- n0-n
for(v in rev_items) x[[v]] <- 7-x[[v]]
stopifnot(!anyNA(x),all(unlist(x)>=1 & unlist(x)<=6)); write.csv(x,"bfi-preprocessing.csv",row.names=FALSE)
pc <- prcomp(x,center=TRUE,scale.=TRUE); eig <- pc$sdev^2; prop <- eig/sum(eig); cum5 <- sum(prop[1:5])
pl <- unclass(varimax(sweep(pc$rotation[,1:5,drop=FALSE],2,pc$sdev[1:5],"*"))$loadings); colnames(pl) <- paste0("RotPC",1:5)
fw <- character(); set.seed(analysis_seed); fa <- withCallingHandlers(factanal(x,factors=5,rotation="varimax",control=list(nstart=10)),warning=function(w){fw<<-c(fw,conditionMessage(w));invokeRestart("muffleWarning")})
fl <- unclass(fa$loadings); colnames(fl) <- paste0("Factor",1:5)
fconv <- if(!is.null(fa$converged)) as.character(fa$converged) else "successful_return_no_explicit_code"; fwarn <- if(length(fw)) paste(unique(fw),collapse=" | ") else "none"
tidy <- function(a,m) data.frame(method=m,item=rep(rownames(a),times=ncol(a)),dimension=rep(colnames(a),each=nrow(a)),loading=as.vector(a))
lo <- rbind(tidy(pl,"PCA"),tidy(fl,"FA")); write.csv(lo,"bfi-loadings.csv",row.names=FALSE)
L <- list(); add <- function(method,metric,dimension="",item="",value=NA_real_,text="",unit="") L[[length(L)+1]] <<- data.frame(method,metric,dimension,item,value,text,unit,stringsAsFactors=FALSE)
add("Data","input_rows",value=n0,unit="persons");add("Data","complete_case_rows",value=n,unit="persons");add("Data","excluded_incomplete_rows",value=nex,unit="persons")
for(v in items)add("Data","missing_count",item=v,value=miss[v],unit="responses")
for(j in 1:5){add("PCA","eigenvalue",paste0("PC",j),value=eig[j]);add("PCA","variance_explained",paste0("PC",j),value=prop[j],unit="proportion");add("PCA","cumulative_variance_explained",paste0("PC",j),value=sum(prop[1:j]),unit="proportion")}
add("PCA","five_component_variance_explained",value=cum5,unit="proportion")
for(v in items)add("FA","uniqueness",item=v,value=unname(fa$uniquenesses[v]),unit="proportion")
add("FA","likelihood_ratio_statistic",value=unname(fa$STATISTIC));add("FA","degrees_of_freedom",value=unname(fa$dof));add("FA","p_value",value=unname(fa$PVAL));add("FA","convergence",text=fconv);add("FA","warnings",text=fwarn);add("FA","random_seed",value=analysis_seed);add("Environment","R_version",text=R.version.string);add("Input","runtime_md5",text=runtime_md5(input_file))
mo <- do.call(rbind,L);write.csv(mo,"bfi-method-metrics.csv",row.names=FALSE,na="")
png("bfi-method-comparison.png",width=2400,height=1800,res=300,bg="white",type="cairo")
layout(matrix(1:3,nrow=1),widths=c(1,1.35,1.35));par(mar=c(7,4.2,3.3,1.2),mgp=c(2.4,.7,0),tcl=-.25,family="sans",las=1)
co <- c("#4477AA","#66AADD","#228833","#CCBB44","#AA3377");bp<-barplot(100*prop[1:5],names.arg=paste0("PC",1:5),col=co,border=NA,ylim=c(0,max(100*prop[1:5])*1.25),ylab="Total variance explained (%)",main="Unrotated PCA variance",cex.names=.75);text(bp,100*prop[1:5],sprintf("%.1f",100*prop[1:5]),pos=3,cex=.78);mtext(sprintf("n = %d complete cases",n),1,4.7,cex=.75)
heat <- function(a,title){zz<-t(a[nrow(a):1,,drop=FALSE]);pal<-colorRampPalette(c("#3B4CC0","#F7F7F7","#B40426"))(101);image(seq_len(ncol(a)),seq_len(nrow(a)),zz,zlim=c(-1,1),col=pal,axes=FALSE,xlab="",ylab="",main=title);axis(1,at=seq_len(ncol(a)),labels=colnames(a),cex.axis=.68,gap.axis=-1);axis(2,at=seq_len(nrow(a)),labels=rev(rownames(a)),cex.axis=.62);box(col="grey55");abline(v=seq(1.5,ncol(a)-.5,1),h=seq(1.5,nrow(a)-.5,1),col="white",lwd=.35)}
heat(pl,"Varimax-rotated PCA loading columns");mtext("Blue negative, red positive; centered at zero",1,4.7,cex=.72);heat(fl,"Varimax ML factor loadings");mtext("Same respondents and item scoring",1,4.7,cex=.72);dev.off()
ih <- runtime_md5(input_file); pcl <- paste(sprintf("- PC%d: %.2f%% (cumulative %.2f%%)",1:5,100*prop[1:5],100*cumsum(prop[1:5])),collapse="\n"); mt <- paste(sprintf("%s=%d",items,miss),collapse=", "); ut <- paste(sprintf("%s=%.3f",items,fa$uniquenesses[items]),collapse=", ")
report <- c("# PCA and exploratory factor analysis of the psych::bfi items","","## Data and reproducibility","",sprintf("- Input: %s",input_file),sprintf("- Runtime input MD5: %s",ih),sprintf("- Download identity SHA-256 (independently verified for the supplied source CSV): %s",source_sha256),sprintf("- Fixed random seed before factanal: %d",analysis_seed),sprintf("- R version: %s",R.version.string),sprintf("- Input rows: %d; complete cases across the 25 items: %d; excluded incomplete rows: %d.",n0,n,nex),sprintf("- Item-level missing counts: %s.",mt),"- Analysis columns: A1-A5, C1-C5, E1-E5, N1-N5, and O1-O5 only. Row identifiers and demographics were not analyzed.","- Reverse scoring used 7 - x for A1, C4, C5, E1, E2, O2, and O5 before either method.","","## PCA results","","PCA used standardized variables. Percentages labeled PC1-PC5 are calculated from the unrotated eigenvalues. The first five eigenvector × square-root-eigenvalue loading columns were then rotated with orthogonal varimax and labeled RotPC1-RotPC5. Rotation redistributes variance among these loading columns but does not change their combined explained variance.","",pcl,"",sprintf("Together, the five retained components explain %.2f%% of total standardized item variance.",100*cum5),"","## Maximum-likelihood factor analysis results","","The five-factor model used maximum-likelihood factanal with orthogonal varimax rotation on the same complete cases.","",sprintf("- Likelihood-ratio statistic: %.3f on %d df; p = %.4g.",unname(fa$STATISTIC),unname(fa$dof),unname(fa$PVAL)),"- Numerical convergence only means the optimizer returned a solution; it does not establish adequate model fit. The very small likelihood-ratio p-value rejects exact five-factor model fit under the maximum-likelihood assumptions.",sprintf("- Optimizer status: %s; warnings: %s.",fconv,fwarn),sprintf("- Item uniquenesses: %s.",ut),"","Uniquenesses estimate item-specific plus error variance outside the five common factors. FA sums of squared loadings are not PCA variance explained: common-factor analysis models shared covariance, whereas PCA components partition total observed variance.","","## Interpretation and limitations","","PCA gives weighted composites summarizing total standardized variance. Maximum-likelihood factor analysis posits latent common factors accounting for covariance. Similar loading patterns can support a shared descriptive interpretation, but do not make the methods equivalent.","","Dimension order and signs are arbitrary after extraction and rotation: a whole component or factor may be multiplied by -1 without changing the solution, and columns can be permuted. Compare loading patterns, not signs or column numbers mechanically.","","The 1-6 responses are ordinal. Pearson correlations, standardized PCA, and normal-theory ML factor analysis treat them as approximately continuous; this pragmatic approximation may differ from ordinal-threshold and polychoric-correlation analyses. Complete-case analysis discards respondents missing any item and can be biased if missingness relates to responses or participant characteristics. This worked example is not a diagnostic instrument or a basis for individual psychological assessment.","","## Output guide","","- bfi-preprocessing.csv: complete-case, reverse-scored 25-item matrix.","- bfi-loadings.csv: tidy rotated loadings for PCA and factor analysis.","- bfi-method-metrics.csv: sample sizes, missingness, unrotated PCA variance quantities, FA uniquenesses and fit/convergence metadata, seed, R version, and runtime input MD5.","- bfi-method-comparison.png: PCA variance and side-by-side loading heatmaps.")
writeLines(report,"bfi-method-report.md",useBytes=TRUE)
files<-c("bfi-method-comparison.R","bfi-preprocessing.csv","bfi-loadings.csv","bfi-method-metrics.csv","bfi-method-comparison.png","bfi-method-report.md");stopifnot(all(file.exists(files)))
pr<-read.csv("bfi-preprocessing.csv",check.names=FALSE);lr<-read.csv("bfi-loadings.csv");mr<-read.csv("bfi-method-metrics.csv");rr<-readLines("bfi-method-report.md",warn=FALSE)
stopifnot(nrow(pr)==n,ncol(pr)==25,nrow(lr)==250,all(table(lr$method)==125),nrow(mr)==77,length(rr)>20,file.info("bfi-method-comparison.png")$size>10000)
cat(sprintf("input=%d complete=%d excluded=%d pca5=%.6f\n",n0,n,nex,cum5));cat(sprintf("fa_stat=%.6f df=%d p=%.8g convergence=%s warnings=%s\n",unname(fa$STATISTIC),unname(fa$dof),unname(fa$PVAL),fconv,fwarn));cat(sprintf("rows_pre=%d rows_load=%d rows_metrics=%d report_lines=%d png_bytes=%d\n",nrow(pr),nrow(lr),nrow(mr),length(rr),file.info("bfi-method-comparison.png")$size));cat(sprintf("runtime_md5=%s\nR=%s\n",ih,R.version.string))