# Programm zur Durchführung einer Faktorenanalyse (2 Faktoren)

# V 0.8 vom 31.10.2006/25.05.2007

# (c) Dr. Alexander Preuß

# Eine vollständige Datenmatrix X muss vorhanden sein!



# Bestimmung von I und J

I <- nrow(X)

J <- ncol(X)

# Bestimmung min(I,J)

r = min(I,J)

# Bestimmung eines Vektors EinsI und EinsJ

EinsI <- c(rep(1, times=I))

EinsJ <- c(rep(1, times=J))

# Bestimmung der Spaltenmittel

smittel <- c(rep(0, times = J))
sstabw <- c(rep(0, times = J))

for (j in 1:J)

	{smittel[j] <- (EinsI%*%X[1:I,j])/I}

# Bestimmung der Spaltenstandardabweichung

for (j in 1:J)

	{

	for (i in 1:I)

		{sstabw[j] <- sstabw[j] + (X[i,j]-smittel[j])^2}

	sstabw[j] <- sqrt(sstabw[j]/(I-1))
	
	}

# Bestimmung der Matrix Z (Standardwerte)

Z <- X

for (i in 1:I)

	{

	for (j in 1:J)

		{Z[i,j] <- (X[i,j]-smittel[j])/sstabw[j]}

	}

# Bestimmung der Korrelationsmatrix R

R <- t(Z)%*%Z/(I-1)



# Bestimmung der Eigenwerte und Eigenvektoren von R

L <- eigen(R)$values
E <- eigen(R)$vectors

# Bestimmung der Eigenwerte > 1 (Kaiser-Kriterium)

Kf <- 0

for (j in 1:J)

	{if (L[j] > 1) Kf = Kf + 1}


# Hier: Automatisch festgelegt auf zwei Faktoren

K <- 2

# Bestimmung der Faktorenladungsmatrix A

A <- matrix(c(rep(0, times = J*K)),nrow=J)

Ka <- c(rep(0, times = J))

A[1:J,1:K] <- E[1:J,1:K]%*%diag(sqrt(L[1:K]))


Rr <- A%*%t(A)

Rest <- R - Rr


A2 <- A[1:J,1:2]

A2t <- t(A2)%*%A2


# Kommunalitäten

Ka = diag(A%*%t(A))	


varsum <- 0

winkel <- 0

# Rotation (Varimax, mit Kaiser-Normalisierung)

for (w in -4500:4500)

	{

	w2 <- w/18000*pi
	
	Rt <- matrix(c(cos(w2), -sin(w2), sin(w2), cos(w2)), byrow=T, nrow=2)

	# Bestimmung der Matrix der quadrierten Ladungen mit Kaiser-Normalisierung

	A2r <- diag(1/sqrt(Ka))%*%A2%*%Rt

	H <- A2r * A2r

	# Bestimmung der Spaltenmittel

	smittelh <- (1/J)*EinsJ%*%H


	# Bestimmung der Spaltenvarianz

		# Abweichungen der Werte von H vom Mittel

		Habw <- H - EinsJ%*%smittelh

		svarh <- Habw*Habw/(J-1)

		varsum2 <- sum(sum(svarh))

		

	if (varsum2 > varsum)

	{

	winkel <- w/100
		
	varsum <- varsum2

	varianz <- sum(sum(svarh))

	mittel <- smittelh

	HH <- H
	
	}				

	}

	
Rt <- matrix(c(cos(winkel/180*pi), -sin(winkel/180*pi), sin(winkel/180*pi), cos(winkel/180*pi)), byrow=T, nrow=2)

Ar <- A2%*%Rt

fak <- c("Faktor 1", "Faktor 2")

Fu <- matrix(c(rep(0,times=J*2)),nrow = J, dimnames = list(cnamen,fak))



for (j in 1:J)

	{for (k in 1:2)

		{

			if (abs(A[j,k]) < 0.5) Fu[j,k] = 0 else Fu[j,k] = A[j,k]

		}
	}



F <- matrix(c(rep(0,times=J*2)),nrow = J, dimnames = list(cnamen,fak))



for (j in 1:J)

	{for (k in 1:2)

		{

			if (abs(Ar[j,k]) < 0.5) F[j,k] = 0 else F[j,k] = Ar[j,k]

		}
	}

# Bestimmung der Anpassungsgüte

Summe_L = sum(L)
G1 = L[1]/Summe_L
G2 = L[2]/Summe_L
G = G1 + G2

# Graphische Ausgabe der Faktorenmatrix

windows()


plot(A[1:J,1],A[1:J,2], pch=".", xlim=c(-1,1), ylim=c(-1,1), xlab="1. Achse", ylab="2. Achse")

abline(h=0)
abline(v=0)
title("Darstellung der Faktorenladungsmatrix")

for (j in 1:J) {text(A[j,1],A[j,2],cnamen[j])}


# Graphische Ausgabe der rotierten Faktorenmatrix

windows()

Kx <- matrix(c(rep(0,times=4*J)),ncol=2)



plot(Ar[1:J,1],Ar[1:J,2], pch=".", xlim=c(-1,1), ylim=c(-1,1), xlab="1. Achse", ylab="2. Achse")


abline(h=0)
abline(v=0)

Rtt = t(Rt)

par(new=T)
plot(c(0,Rtt[1,1]*1000), c(0,Rtt[2,1]*1000),type="l", xlim=c(-1,1), ylim=c(-1,1), col="grey", xlab = "", ylab = "")

par(new=T)
plot(c(0,-Rtt[1,1]*1000), c(0,-Rtt[2,1]*1000),type="l", xlim=c(-1,1), ylim=c(-1,1), col="grey", xlab = "", ylab = "")


par(new=T)
plot(c(0,Rtt[1,2]*1000), c(0,Rtt[2,2]*1000),type="l", xlim=c(-1,1), ylim=c(-1,1), col="grey", xlab = "", ylab = "")

par(new=T)
plot(c(0,-Rtt[1,2]*1000), c(0,-Rtt[2,2]*1000),type="l", xlim=c(-1,1), ylim=c(-1,1), col="grey", xlab = "", ylab = "")



title("Darstellung der rotierten Faktorenladungsmatrix,
in grau: Achsen des alten Koordinatensystems")

for (j in 1:J) {text(Ar[j,1],Ar[j,2],cnamen[j])}


# Screeplot

zaehler <- c(rep(0, times=J))

for (j in 1:J)

	{zaehler[j] = j}

windows()

plot(zaehler[1:J],L[1:J], pch="*", xlab = "Faktor", ylab = "Eigenwert")
title("Screeplot")
abline(h=1)


print (" ")
print (" ")
print (" ")
print ("Ergebnis der Faktorenanalyse")
print (" ")
print (" ")
print("Korrelationsmatrix")
print (" ")
print(R)
print (" ")
print (" ")
print("Zahl der Faktoren (Kaiser-Kriterium)")
print (" ")
print(Kf)
print (" ")
print (" ")
print("Kommunalitäten (nach Extraktion)")
print (" ")
print (Ka)
print (" ")
print (" ")
print ("Komponentenmatrix")
print (" ")
print (A)
print (" ")
print (" ")
print ("rotierte Komponentenmatrix (bei zwei Faktoren)")
print (" ")
print (Ar)
print (" ")
print (" ")
print ("rotierte Komponentenmatrix, ohne Werte < 0,5")
print (" ")
print (F)
print (" ")
print (" ")
print ("Rotationsmatrix")
print (" ")
print (Rt)
print (" ")
print (" ")
print ("Rotationswinkel")
print (" ")
print (winkel)