0.1 Iran

0.2 Describe the aim and the research question of your small study - they might be different from those described in the article (1 pt)

This study is part of a huge imaginary project which aimed at testing TIMSS for different countries.

“How universal are these scales and what are the differences?”

• Test if western created measures works in eastern Arab country
• Test explanatory power of attitudes factors to math achievements

2. Describe the data, give descriptive statistics for your variables - should be the same as in the study (2 pts)

library(foreign)

BSGIRNM6 <- read.spss("~/datanal/3year/FA/BSGIRNM6.sav", to.data.frame=TRUE)

library(dplyr)
df_Iran <- BSGIRNM6 %>% dplyr::select(BSBM19A:BSBM19D,  BSBM17A, BSBM20A:BSBM20E, BSMMAT01, ITSEX, BSBG10A, BSBG07B, BSBG07A)

knitr::kable(head(df_Iran, 5))

code.book <- data.frame(code = colnames(df_Iran), question = c(" I usually do well in mathematics","Mathematics is more difficult for me than for many of my classmates",
"Mathematics is not one of my strengths", " I learn things quickly in mathematics", " I enjoy learning mathematics", " I think learning mathematics will help me in my daily life",
" I need mathematics to learn other school subjects", "I need to do well in mathematics to get into the <university> of my choice", " I need to do well in mathematics to get the job I want", " I would like a job that involves using mathematics", " the eighth grade mathematics achievement plausible values", "sex", "Were you born in Iran?", "What is the highest level of education completed by your father", "What is the highest level of education completed by your mother"))
knitr::kable(code.book)

library(psych)
library(formattable)

formattable(describe(df_Iran))

Mean age is near 18

a <-  nrow(df_Iran)
b <-  nrow(na.omit(df_Iran))
cat("There were", a, "observations in collected by TIMSS data.", "Without missings that would be", b)

## There were 6130 observations in collected by TIMSS data. Without missings that would be 5580

library(ggplot2)
library(ggthemes)
library(gridExtra)
df_Iran_o <- df_Iran %>% na.omit %>% select(-BSMMAT01)
  
df_Iran_reg <- df_Iran %>% na.omit %>% select(BSMMAT01)
plots <- list()  

for (nm in names(df_Iran_o)){
 plots[[nm]] <-  ggplot(df_Iran_o)+
   ggtitle(paste("distribution of variable" ,nm))+
    geom_bar(aes_string(nm), stat = "count", fill = "#E63946")+ 
    theme_wsj() + 
    coord_flip()+
   theme(plot.title = element_text(size = 18, face = "bold"))
}


plots[[sample(1:12, 1)]]

table(df_Iran_o$ITSEX)

## 
## Female   Male 
##   2752   2828

There are 2752 female and 2828 male students in the sample.

0.3 Conduct EFA. Present your correlation matrix. How did you choose the number of factors? (1 pt)

library(polycor)

corrr<- hetcor(df_Iran_o[,-10:-14])
corrr <- corrr$correlations


cor.mtest <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], ...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}
# matrix of the p-value of the correlation
p.mat <- cor.mtest(corrr)


corrplot::corrplot(corrr, type="upper", order="hclust", 
         p.mat = p.mat, sig.level = 0.05, insig = "blank")

All variables shows significant p-value < 0.05 correlation Most of them are strongly inter-correlated (almost 0.7) There are two group of correlated variables

fa.parallel(corrr, fa="both", n.iter=100, n.obs = 5580) 

## Parallel analysis suggests that the number of factors =  4  and the number of components =  2

4 factors were offered, so lets try them

fa.diagram(fa(corrr, nfactors=4, rotate="varimax", fm="ml", n.obs = 5577))

There are two factors, that explains only 2 variables each. We need more explanation!

fa.diagram(fa(corrr, nfactors=3, rotate="varimax", fm="ml", n.obs = 5577))

Still only two for the third factor.

fa.diagram(fa(corrr, nfactors=2, rotate="varimax", fm="ml", n.obs = 5577))

Better with 2 factors, even if one loading is not good. (0.4)

0.4 Interpret the resulting factor structure. What type of rotation did you choose and why? Show the diagram. Is your structure different from the one Hanna Eklof got in her article? Why? (2 pts)

Varimax was used as rotation, because oblimin showed correlation between factor, when other variables does not changed much. rotation was used to significantly increase quality of model, to relay on data not on random. Low correlation 0.4 was stable with oblimin rotation to, so I decided to not follow the article.

BSBM19B, BSBM19C, BSBM17A, BSBM19D, BSBM19A – ML2

BSBM20A, BSBM20B, BSBM20D, BSBM20C – ML1

ML1 match is necessary for human in future

ML2 skills in match (19c and 19b correspond lack of skills)

Factors do meet results of PCA from the paper.

1 Analyze the model fit; Chronbach’s alpha. What can you say about your scale(s)? (1 pt)

code.book$code

##  [1] BSBM19A  BSBM19B  BSBM19C  BSBM19D  BSBM17A  BSBM20A  BSBM20B 
##  [8] BSBM20C  BSBM20D  BSBM20E  BSMMAT01 ITSEX    BSBG10A  BSBG07B 
## [15] BSBG07A 
## 15 Levels: BSBG07A BSBG07B BSBG10A BSBM17A BSBM19A BSBM19B ... ITSEX

df_Iran_num <- sapply( df_Iran_o, as.numeric )
df_Iran_num <- as.data.frame(df_Iran_num)
Iran_fa_model <- fa(df_Iran_num[,-11:-14], nfactors=2, rotate="varimax", fm="ml", n.obs = 5577, cor = "mixed") 

## 
## mixed.cor is deprecated, please use mixedCor.

print(Iran_fa_model$loadings,cutoff = 0.3)

## 
## Loadings:
##         ML1    ML2   
## BSBM19A  0.354  0.774
## BSBM19B        -0.397
## BSBM19C        -0.503
## BSBM19D  0.338  0.744
## BSBM17A  0.451  0.652
## BSBM20A  0.571  0.362
## BSBM20B  0.571       
## BSBM20C  0.880       
## BSBM20D  0.896       
## BSBM20E  0.646  0.405
## 
##                  ML1   ML2
## SS loadings    3.106 2.350
## Proportion Var 0.311 0.235
## Cumulative Var 0.311 0.546

BSBM20E and BSBM17A quite fiited in both factors

summary(Iran_fa_model) 

## 
## Factor analysis with Call: fa(r = df_Iran_num[, -11:-14], nfactors = 2, n.obs = 5577, rotate = "varimax", 
##     fm = "ml", cor = "mixed")
## 
## Test of the hypothesis that 2 factors are sufficient.
## The degrees of freedom for the model is 26  and the objective function was  0.53 
## The number of observations was  5580  with Chi Square =  2950.72  with prob <  0 
## 
## The root mean square of the residuals (RMSA) is  0.07 
## The df corrected root mean square of the residuals is  0.09 
## 
## Tucker Lewis Index of factoring reliability =  0.828
## RMSEA index =  0.142  and the 10 % confidence intervals are  0.138 0.146
## BIC =  2726.42

RMSA is bigger then 0.5, same as RMSR more then 0.05, both of them does not show well model fit. Tucker Lewis Index does not cross 0.9 border. Not a good model

psych::alpha(df_Iran_num[c('BSBM19A',  'BSBM19B', 'BSBM19C',  'BSBM19D',  'BSBM17A',  'BSBM20A',  'BSBM20B',  'BSBM20C',  'BSBM20D', 'BSBM20E')], check.keys = TRUE)

## 
## Reliability analysis   
## Call: psych::alpha(x = df_Iran_num[c("BSBM19A", "BSBM19B", "BSBM19C", 
##     "BSBM19D", "BSBM17A", "BSBM20A", "BSBM20B", "BSBM20C", "BSBM20D", 
##     "BSBM20E")], check.keys = TRUE)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd median_r
##       0.82      0.83    0.85      0.32 4.7 0.0036  1.9 0.59     0.34
## 
##  lower alpha upper     95% confidence boundaries
## 0.81 0.82 0.83 
## 
##  Reliability if an item is dropped:
##          raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## BSBM19A       0.79      0.80    0.82      0.30 3.9   0.0043 0.030  0.34
## BSBM19B-      0.84      0.84    0.85      0.37 5.3   0.0033 0.018  0.36
## BSBM19C-      0.82      0.82    0.84      0.34 4.7   0.0036 0.029  0.36
## BSBM19D       0.79      0.80    0.83      0.31 4.0   0.0042 0.030  0.34
## BSBM17A       0.79      0.80    0.83      0.30 3.9   0.0043 0.030  0.31
## BSBM20A       0.80      0.81    0.83      0.32 4.2   0.0040 0.030  0.31
## BSBM20B       0.81      0.81    0.84      0.33 4.4   0.0039 0.030  0.34
## BSBM20C       0.80      0.81    0.83      0.32 4.3   0.0039 0.027  0.34
## BSBM20D       0.80      0.81    0.82      0.32 4.2   0.0040 0.027  0.34
## BSBM20E       0.79      0.80    0.83      0.31 4.0   0.0043 0.030  0.31
## 
##  Item statistics 
##             n raw.r std.r r.cor r.drop mean   sd
## BSBM19A  5580  0.72  0.72  0.70   0.64  2.0 0.92
## BSBM19B- 5580  0.39  0.35  0.25   0.22  2.5 1.07
## BSBM19C- 5580  0.53  0.50  0.42   0.38  2.5 1.11
## BSBM19D  5580  0.70  0.70  0.67   0.60  2.0 0.95
## BSBM17A  5580  0.73  0.73  0.69   0.63  1.8 0.97
## BSBM20A  5580  0.63  0.65  0.60   0.53  1.6 0.85
## BSBM20B  5580  0.58  0.60  0.53   0.47  1.8 0.88
## BSBM20C  5580  0.60  0.63  0.59   0.50  1.5 0.81
## BSBM20D  5580  0.63  0.65  0.63   0.53  1.6 0.88
## BSBM20E  5580  0.72  0.71  0.68   0.61  2.3 1.09
## 
## Non missing response frequency for each item
##            1    2    3    4 miss
## BSBM19A 0.33 0.42 0.16 0.09    0
## BSBM19B 0.19 0.35 0.21 0.26    0
## BSBM19C 0.21 0.32 0.19 0.28    0
## BSBM19D 0.34 0.38 0.18 0.10    0
## BSBM17A 0.49 0.31 0.10 0.09    0
## BSBM20A 0.61 0.26 0.07 0.06    0
## BSBM20B 0.48 0.34 0.11 0.06    0
## BSBM20C 0.69 0.21 0.06 0.05    0
## BSBM20D 0.60 0.24 0.09 0.06    0
## BSBM20E 0.31 0.32 0.19 0.19    0

Cronbach’s Alpha shows good inter-item reliability, because value is near 0.82, (expected 0.9 or more, for exellent) SD is nice but, who cares.

load <- Iran_fa_model$loadings[,1:2] 
plot(load) # set up plot 

fascores<-as.data.frame(Iran_fa_model$scores)
Iran_fa<-cbind(df_Iran_o,fascores)
Iran_fa$BSMMAT01 <- df_Iran_reg$BSMMAT01

model1 <- lm(as.numeric(BSMMAT01) ~ ML1 + ML2 + ITSEX + BSBG10A + ordered(BSBG07B) + ordered(BSBG07A), data = Iran_fa)
library("QuantPsyc")

summary(model1)

## 
## Call:
## lm(formula = as.numeric(BSMMAT01) ~ ML1 + ML2 + ITSEX + BSBG10A + 
##     ordered(BSBG07B) + ordered(BSBG07A), data = Iran_fa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4764.4 -1117.6    80.9  1081.7  4562.9 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3422.692     32.737 104.552  < 2e-16 ***
## ML1                 -28.331     19.683  -1.439 0.150108    
## ML2                -656.846     21.532 -30.505  < 2e-16 ***
## ITSEXMale          -144.673     39.169  -3.694 0.000223 ***
## BSBG10ANo          -702.698    208.502  -3.370 0.000756 ***
## ordered(BSBG07B).L  619.213     82.667   7.490 7.94e-14 ***
## ordered(BSBG07B).Q -660.461     78.051  -8.462  < 2e-16 ***
## ordered(BSBG07B).C -588.181     73.229  -8.032 1.16e-15 ***
## ordered(BSBG07B)^4 -320.505     69.324  -4.623 3.86e-06 ***
## ordered(BSBG07B)^5 -352.518     64.049  -5.504 3.88e-08 ***
## ordered(BSBG07B)^6 -151.424     59.202  -2.558 0.010562 *  
## ordered(BSBG07B)^7  192.405     59.661   3.225 0.001267 ** 
## ordered(BSBG07A).L  365.225     89.491   4.081 4.54e-05 ***
## ordered(BSBG07A).Q -841.559     78.946 -10.660  < 2e-16 ***
## ordered(BSBG07A).C -489.827     80.254  -6.103 1.11e-09 ***
## ordered(BSBG07A)^4 -251.322     82.666  -3.040 0.002375 ** 
## ordered(BSBG07A)^5 -224.154     76.689  -2.923 0.003482 ** 
## ordered(BSBG07A)^6    4.865     66.724   0.073 0.941881    
## ordered(BSBG07A)^7  288.831     61.575   4.691 2.79e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1449 on 5561 degrees of freedom
## Multiple R-squared:  0.3263, Adjusted R-squared:  0.3241 
## F-statistic: 149.6 on 18 and 5561 DF,  p-value: < 2.2e-16

Iran_fa_num <- sapply( Iran_fa, as.numeric )
Iran_fa_num <- as.data.frame(
  Iran_fa_num
)

model2 <- lm(as.numeric(BSMMAT01) ~ ML1 + ML2 + ITSEX + BSBG10A + BSBG07B + BSBG07A, data = Iran_fa_num)
library("QuantPsyc")

lm.beta(model2)

##          ML1          ML2        ITSEX      BSBG10A      BSBG07B 
## -0.005663678 -0.383103322 -0.054519928 -0.044915069  0.167236755 
##      BSBG07A 
##  0.130957784