This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

rm(list = ls())

#loading data
d <-read.csv(file.choose(), header = T)
head(d)

#lets look at the structure
str(d)
#summary
summary(d)

#Problem statment states respondents used a 5 point likert scale
#there are few variables which has values as 6
#We will convert them as 5 as 5 being the highest factor

d[d==6] <- 5  
d[,c(12,25)] <- 6 - d[,c(12,25)]
dim(d)

library(psych)
library(corrplot)
library(nFactors)
library(MASS)
library(GPArotation)
#calling libraries


corfac <- cor(d[-1])
corrplot(corfac)

#To do PCA or Factor; first lets do the KMO and Bartlett test

#preforming KMO test
KMO(d[,-1])

#overall MSA is 0.85 which is good enough for the model
#FACTOR analysis can be carried out for this data

#performing barllet test

print(cortest.bartlett(corfac, nrow(d)))
#P value is significantly low and hence we can assume
#there is some significiance b/w the variables


#lets do a PCA to identify no of factors to be used


pca1 <- princomp(d[,-1], scores = T, cor = T)
pca1
summary(pca1)
#from the summary, it looks like there are 4 components
#SD of upto comp 5 is greater than 1 (eigen values)
#It explains upto 62% of the model

screeplot(pca1, type = "b", main = "Scree Plot of PCA")

pca2 <- principal(d[,-1], nfactors = length(d[,-1]), rotate = "none")
pca2
cat("\n\n Eigen values \n")
print(pca2$values)
#Even in this principal function 5 components are significant
plot(pca2$values, type = "b",xlab = "Factors", ylab = "Eigen values", main = "SCREE PLOT") 

#alternative
no <- fa.parallel(d[,-1])


#We will use 5 factors for our factor analysis

sol1 <- fa(r=d[,-1], nfactors = 5, rotate = "oblimin", fm="pa")
print(sol1)
sol1$uniquenesses

#factor 5 explains only 0.1 variance
#lets try using 4 factors
#using oblimin rotation as there is some correlation within the factors as shown in corrplot
#varimax & promax also can be used
sol2 <- fa(d[,-1], nfactors = 4, rotate = "oblimin", fm="ml")
print(sol2)

#the loadings are distributed
#using a cutoff value of 0.3, lets see the factors


print(sol2$loadings, cutoff = 0.3)
#filling, natural, fiber & health are contributing to factor 1
#sugar, salt & sweet are in factor 2
#kids, family & economical are in factor 3
#plain has a negative correlation in factor 4


#Factor loading can be classified based on their magnitude.
#Value  Interpretation
#> 30   Minimum consideration level
#> 40   More important
#> 50   Practically significant

fa.diagram(sol2)
#This diagram tell us the pictures

fl <- round(unclass(sol2$loadings),2)
fl
write.csv(fl, "Factors for ds.csv")


library(lavaan)
#4 factors adding them to one unique variable 
factor1 <- c(2,3,4,8,9,14,19,23,26)
factor2 <- c(5,7,16,20,22)
factor3 <- c(11,13,15)
factor4 <- c(10,12,17,18,21,24,25)

models <- list()
fits <- list()


models$m3 <- 
  ' Health =~ Filling + Natural + Fibre + Satisfying +Energy +Regular +Quality +Nutritious
    Taste =~ Sweet + Salt + Calories + Sugar +Process 
    Famaly =~ Kids + Economical + Family 
    Result =~ Fun + Soggy + Plain + Crisp + Fruit + Treat + Boring 
  '


fits$m3 <- lavaan::cfa(models$m3, data = d)
fits$m3
standardizedSolution(fits$m3)

summary(fits$m3)





#Crombach's alpha is a measure of internal consistency, 
#that is, how closely related a set of items are as a group.


factor1alpha <- psych::alpha(d[,factor1], check.keys = TRUE)
factor2alpha <- psych::alpha(d[,factor2], check.keys = TRUE)
factor3alpha <- psych::alpha(d[,factor3], check.keys = TRUE)
factor4alpha <- psych::alpha(d[,factor4], check.keys = TRUE)


factor1alpha$total$raw_alpha
factor2alpha$total$raw_alpha
factor3alpha$total$raw_alpha
factor4alpha$total$raw_alpha

#alpha values are 0.7

d[,factor1]
#Just looking at factor 1 if it is alinged

#better to take the average of the factors and assign them to factors

?apply
#tried with table function; but mean is better 

#lets add all the factors to the original data set

d$factor1Score <- apply(d[,factor1],1,mean)
d$factor2Score <- apply(d[,factor2],1,mean)
d$factor3Score <- apply(d[,factor3],1,mean)
d$factor4Score <- apply(d[,factor4],1,mean)
colnames(d)[27:30] <-c("Health", "Taste", "Family", "Result")
?aggregate


#lets take our the delimited factors
d_full<-aggregate(d[,27:30],  list(d[,1]), mean)



#final factors
format(d_full, digits = 2)
?tapply
tapply(d[,1],d[,27:30],mean)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpUaGlzIGlzIGFuIFtSIE1hcmtkb3duXShodHRwOi8vcm1hcmtkb3duLnJzdHVkaW8uY29tKSBOb3RlYm9vay4gV2hlbiB5b3UgZXhlY3V0ZSBjb2RlIHdpdGhpbiB0aGUgbm90ZWJvb2ssIHRoZSByZXN1bHRzIGFwcGVhciBiZW5lYXRoIHRoZSBjb2RlLiANCg0KVHJ5IGV4ZWN1dGluZyB0aGlzIGNodW5rIGJ5IGNsaWNraW5nIHRoZSAqUnVuKiBidXR0b24gd2l0aGluIHRoZSBjaHVuayBvciBieSBwbGFjaW5nIHlvdXIgY3Vyc29yIGluc2lkZSBpdCBhbmQgcHJlc3NpbmcgKkN0cmwrU2hpZnQrRW50ZXIqLiANCg0KYGBge3J9DQpybShsaXN0ID0gbHMoKSkNCg0KI2xvYWRpbmcgZGF0YQ0KZCA8LXJlYWQuY3N2KGZpbGUuY2hvb3NlKCksIGhlYWRlciA9IFQpDQpoZWFkKGQpDQoNCiNsZXRzIGxvb2sgYXQgdGhlIHN0cnVjdHVyZQ0Kc3RyKGQpDQojc3VtbWFyeQ0Kc3VtbWFyeShkKQ0KDQojUHJvYmxlbSBzdGF0bWVudCBzdGF0ZXMgcmVzcG9uZGVudHMgdXNlZCBhIDUgcG9pbnQgbGlrZXJ0IHNjYWxlDQojdGhlcmUgYXJlIGZldyB2YXJpYWJsZXMgd2hpY2ggaGFzIHZhbHVlcyBhcyA2DQojV2Ugd2lsbCBjb252ZXJ0IHRoZW0gYXMgNSBhcyA1IGJlaW5nIHRoZSBoaWdoZXN0IGZhY3Rvcg0KDQpkW2Q9PTZdIDwtIDUgIA0KZFssYygxMiwyNSldIDwtIDYgLSBkWyxjKDEyLDI1KV0NCmRpbShkKQ0KDQpsaWJyYXJ5KHBzeWNoKQ0KbGlicmFyeShjb3JycGxvdCkNCmxpYnJhcnkobkZhY3RvcnMpDQpsaWJyYXJ5KE1BU1MpDQpsaWJyYXJ5KEdQQXJvdGF0aW9uKQ0KI2NhbGxpbmcgbGlicmFyaWVzDQoNCg0KY29yZmFjIDwtIGNvcihkWy0xXSkNCmNvcnJwbG90KGNvcmZhYykNCg0KI1RvIGRvIFBDQSBvciBGYWN0b3I7IGZpcnN0IGxldHMgZG8gdGhlIEtNTyBhbmQgQmFydGxldHQgdGVzdA0KDQojcHJlZm9ybWluZyBLTU8gdGVzdA0KS01PKGRbLC0xXSkNCg0KI292ZXJhbGwgTVNBIGlzIDAuODUgd2hpY2ggaXMgZ29vZCBlbm91Z2ggZm9yIHRoZSBtb2RlbA0KI0ZBQ1RPUiBhbmFseXNpcyBjYW4gYmUgY2FycmllZCBvdXQgZm9yIHRoaXMgZGF0YQ0KDQojcGVyZm9ybWluZyBiYXJsbGV0IHRlc3QNCg0KcHJpbnQoY29ydGVzdC5iYXJ0bGV0dChjb3JmYWMsIG5yb3coZCkpKQ0KI1AgdmFsdWUgaXMgc2lnbmlmaWNhbnRseSBsb3cgYW5kIGhlbmNlIHdlIGNhbiBhc3N1bWUNCiN0aGVyZSBpcyBzb21lIHNpZ25pZmljaWFuY2UgYi93IHRoZSB2YXJpYWJsZXMNCg0KDQojbGV0cyBkbyBhIFBDQSB0byBpZGVudGlmeSBubyBvZiBmYWN0b3JzIHRvIGJlIHVzZWQNCg0KDQpwY2ExIDwtIHByaW5jb21wKGRbLC0xXSwgc2NvcmVzID0gVCwgY29yID0gVCkNCnBjYTENCnN1bW1hcnkocGNhMSkNCiNmcm9tIHRoZSBzdW1tYXJ5LCBpdCBsb29rcyBsaWtlIHRoZXJlIGFyZSA0IGNvbXBvbmVudHMNCiNTRCBvZiB1cHRvIGNvbXAgNSBpcyBncmVhdGVyIHRoYW4gMSAoZWlnZW4gdmFsdWVzKQ0KI0l0IGV4cGxhaW5zIHVwdG8gNjIlIG9mIHRoZSBtb2RlbA0KDQpzY3JlZXBsb3QocGNhMSwgdHlwZSA9ICJiIiwgbWFpbiA9ICJTY3JlZSBQbG90IG9mIFBDQSIpDQoNCnBjYTIgPC0gcHJpbmNpcGFsKGRbLC0xXSwgbmZhY3RvcnMgPSBsZW5ndGgoZFssLTFdKSwgcm90YXRlID0gIm5vbmUiKQ0KcGNhMg0KY2F0KCJcblxuIEVpZ2VuIHZhbHVlcyBcbiIpDQpwcmludChwY2EyJHZhbHVlcykNCiNFdmVuIGluIHRoaXMgcHJpbmNpcGFsIGZ1bmN0aW9uIDUgY29tcG9uZW50cyBhcmUgc2lnbmlmaWNhbnQNCnBsb3QocGNhMiR2YWx1ZXMsIHR5cGUgPSAiYiIseGxhYiA9ICJGYWN0b3JzIiwgeWxhYiA9ICJFaWdlbiB2YWx1ZXMiLCBtYWluID0gIlNDUkVFIFBMT1QiKSANCg0KI2FsdGVybmF0aXZlDQpubyA8LSBmYS5wYXJhbGxlbChkWywtMV0pDQoNCg0KI1dlIHdpbGwgdXNlIDUgZmFjdG9ycyBmb3Igb3VyIGZhY3RvciBhbmFseXNpcw0KDQpzb2wxIDwtIGZhKHI9ZFssLTFdLCBuZmFjdG9ycyA9IDUsIHJvdGF0ZSA9ICJvYmxpbWluIiwgZm09InBhIikNCnByaW50KHNvbDEpDQpzb2wxJHVuaXF1ZW5lc3Nlcw0KDQojZmFjdG9yIDUgZXhwbGFpbnMgb25seSAwLjEgdmFyaWFuY2UNCiNsZXRzIHRyeSB1c2luZyA0IGZhY3RvcnMNCiN1c2luZyBvYmxpbWluIHJvdGF0aW9uIGFzIHRoZXJlIGlzIHNvbWUgY29ycmVsYXRpb24gd2l0aGluIHRoZSBmYWN0b3JzIGFzIHNob3duIGluIGNvcnJwbG90DQojdmFyaW1heCAmIHByb21heCBhbHNvIGNhbiBiZSB1c2VkDQpzb2wyIDwtIGZhKGRbLC0xXSwgbmZhY3RvcnMgPSA0LCByb3RhdGUgPSAib2JsaW1pbiIsIGZtPSJtbCIpDQpwcmludChzb2wyKQ0KDQojdGhlIGxvYWRpbmdzIGFyZSBkaXN0cmlidXRlZA0KI3VzaW5nIGEgY3V0b2ZmIHZhbHVlIG9mIDAuMywgbGV0cyBzZWUgdGhlIGZhY3RvcnMNCg0KDQpwcmludChzb2wyJGxvYWRpbmdzLCBjdXRvZmYgPSAwLjMpDQojZmlsbGluZywgbmF0dXJhbCwgZmliZXIgJiBoZWFsdGggYXJlIGNvbnRyaWJ1dGluZyB0byBmYWN0b3IgMQ0KI3N1Z2FyLCBzYWx0ICYgc3dlZXQgYXJlIGluIGZhY3RvciAyDQoja2lkcywgZmFtaWx5ICYgZWNvbm9taWNhbCBhcmUgaW4gZmFjdG9yIDMNCiNwbGFpbiBoYXMgYSBuZWdhdGl2ZSBjb3JyZWxhdGlvbiBpbiBmYWN0b3IgNA0KDQoNCiNGYWN0b3IgbG9hZGluZyBjYW4gYmUgY2xhc3NpZmllZCBiYXNlZCBvbiB0aGVpciBtYWduaXR1ZGUuDQojVmFsdWUJSW50ZXJwcmV0YXRpb24NCiM+IDMwCU1pbmltdW0gY29uc2lkZXJhdGlvbiBsZXZlbA0KIz4gNDAJTW9yZSBpbXBvcnRhbnQNCiM+IDUwCVByYWN0aWNhbGx5IHNpZ25pZmljYW50DQoNCmZhLmRpYWdyYW0oc29sMikNCiNUaGlzIGRpYWdyYW0gdGVsbCB1cyB0aGUgcGljdHVyZXMNCg0KZmwgPC0gcm91bmQodW5jbGFzcyhzb2wyJGxvYWRpbmdzKSwyKQ0KZmwNCndyaXRlLmNzdihmbCwgIkZhY3RvcnMgZm9yIGRzLmNzdiIpDQoNCg0KbGlicmFyeShsYXZhYW4pDQojNCBmYWN0b3JzIGFkZGluZyB0aGVtIHRvIG9uZSB1bmlxdWUgdmFyaWFibGUgDQpmYWN0b3IxIDwtIGMoMiwzLDQsOCw5LDE0LDE5LDIzLDI2KQ0KZmFjdG9yMiA8LSBjKDUsNywxNiwyMCwyMikNCmZhY3RvcjMgPC0gYygxMSwxMywxNSkNCmZhY3RvcjQgPC0gYygxMCwxMiwxNywxOCwyMSwyNCwyNSkNCg0KbW9kZWxzIDwtIGxpc3QoKQ0KZml0cyA8LSBsaXN0KCkNCg0KDQptb2RlbHMkbTMgPC0gDQogICcgSGVhbHRoID1+IEZpbGxpbmcgKyBOYXR1cmFsICsgRmlicmUgKyBTYXRpc2Z5aW5nICtFbmVyZ3kgK1JlZ3VsYXIgK1F1YWxpdHkgK051dHJpdGlvdXMNCiAgICBUYXN0ZSA9fiBTd2VldCArIFNhbHQgKyBDYWxvcmllcyArIFN1Z2FyICtQcm9jZXNzIA0KICAgIEZhbWFseSA9fiBLaWRzICsgRWNvbm9taWNhbCArIEZhbWlseSANCiAgICBSZXN1bHQgPX4gRnVuICsgU29nZ3kgKyBQbGFpbiArIENyaXNwICsgRnJ1aXQgKyBUcmVhdCArIEJvcmluZyANCiAgJw0KDQoNCmZpdHMkbTMgPC0gbGF2YWFuOjpjZmEobW9kZWxzJG0zLCBkYXRhID0gZCkNCmZpdHMkbTMNCnN0YW5kYXJkaXplZFNvbHV0aW9uKGZpdHMkbTMpDQoNCnN1bW1hcnkoZml0cyRtMykNCg0KDQoNCg0KDQojQ3JvbWJhY2gncyBhbHBoYSBpcyBhIG1lYXN1cmUgb2YgaW50ZXJuYWwgY29uc2lzdGVuY3ksIA0KI3RoYXQgaXMsIGhvdyBjbG9zZWx5IHJlbGF0ZWQgYSBzZXQgb2YgaXRlbXMgYXJlIGFzIGEgZ3JvdXAuDQoNCg0KZmFjdG9yMWFscGhhIDwtIHBzeWNoOjphbHBoYShkWyxmYWN0b3IxXSwgY2hlY2sua2V5cyA9IFRSVUUpDQpmYWN0b3IyYWxwaGEgPC0gcHN5Y2g6OmFscGhhKGRbLGZhY3RvcjJdLCBjaGVjay5rZXlzID0gVFJVRSkNCmZhY3RvcjNhbHBoYSA8LSBwc3ljaDo6YWxwaGEoZFssZmFjdG9yM10sIGNoZWNrLmtleXMgPSBUUlVFKQ0KZmFjdG9yNGFscGhhIDwtIHBzeWNoOjphbHBoYShkWyxmYWN0b3I0XSwgY2hlY2sua2V5cyA9IFRSVUUpDQoNCg0KZmFjdG9yMWFscGhhJHRvdGFsJHJhd19hbHBoYQ0KZmFjdG9yMmFscGhhJHRvdGFsJHJhd19hbHBoYQ0KZmFjdG9yM2FscGhhJHRvdGFsJHJhd19hbHBoYQ0KZmFjdG9yNGFscGhhJHRvdGFsJHJhd19hbHBoYQ0KDQojYWxwaGEgdmFsdWVzIGFyZSAwLjcNCg0KZFssZmFjdG9yMV0NCiNKdXN0IGxvb2tpbmcgYXQgZmFjdG9yIDEgaWYgaXQgaXMgYWxpbmdlZA0KDQojYmV0dGVyIHRvIHRha2UgdGhlIGF2ZXJhZ2Ugb2YgdGhlIGZhY3RvcnMgYW5kIGFzc2lnbiB0aGVtIHRvIGZhY3RvcnMNCg0KP2FwcGx5DQojdHJpZWQgd2l0aCB0YWJsZSBmdW5jdGlvbjsgYnV0IG1lYW4gaXMgYmV0dGVyIA0KDQojbGV0cyBhZGQgYWxsIHRoZSBmYWN0b3JzIHRvIHRoZSBvcmlnaW5hbCBkYXRhIHNldA0KDQpkJGZhY3RvcjFTY29yZSA8LSBhcHBseShkWyxmYWN0b3IxXSwxLG1lYW4pDQpkJGZhY3RvcjJTY29yZSA8LSBhcHBseShkWyxmYWN0b3IyXSwxLG1lYW4pDQpkJGZhY3RvcjNTY29yZSA8LSBhcHBseShkWyxmYWN0b3IzXSwxLG1lYW4pDQpkJGZhY3RvcjRTY29yZSA8LSBhcHBseShkWyxmYWN0b3I0XSwxLG1lYW4pDQpjb2xuYW1lcyhkKVsyNzozMF0gPC1jKCJIZWFsdGgiLCAiVGFzdGUiLCAiRmFtaWx5IiwgIlJlc3VsdCIpDQo/YWdncmVnYXRlDQoNCg0KI2xldHMgdGFrZSBvdXIgdGhlIGRlbGltaXRlZCBmYWN0b3JzDQpkX2Z1bGw8LWFnZ3JlZ2F0ZShkWywyNzozMF0sICBsaXN0KGRbLDFdKSwgbWVhbikNCg0KDQoNCiNmaW5hbCBmYWN0b3JzDQpmb3JtYXQoZF9mdWxsLCBkaWdpdHMgPSAyKQ0KP3RhcHBseQ0KdGFwcGx5KGRbLDFdLGRbLDI3OjMwXSxtZWFuKQ0KYGBgDQoNCkFkZCBhIG5ldyBjaHVuayBieSBjbGlja2luZyB0aGUgKkluc2VydCBDaHVuayogYnV0dG9uIG9uIHRoZSB0b29sYmFyIG9yIGJ5IHByZXNzaW5nICpDdHJsK0FsdCtJKi4NCg0KV2hlbiB5b3Ugc2F2ZSB0aGUgbm90ZWJvb2ssIGFuIEhUTUwgZmlsZSBjb250YWluaW5nIHRoZSBjb2RlIGFuZCBvdXRwdXQgd2lsbCBiZSBzYXZlZCBhbG9uZ3NpZGUgaXQgKGNsaWNrIHRoZSAqUHJldmlldyogYnV0dG9uIG9yIHByZXNzICpDdHJsK1NoaWZ0K0sqIHRvIHByZXZpZXcgdGhlIEhUTUwgZmlsZSkuDQo=