library(wordbankr) # WB data
library(tidyverse) # tidy
library(mirt) # IRT models
library(ltm) # more IRT functions
library(psych) # some psychometric stuff (tests of dimensionality)
library(Gifi)# some more psychometric stuff (tests of dimensionality)
library(knitr) # some formatting, tables, etc
library(patchwork) # combining plots. 
library(GGally) # More plottinng options. 
library(lordif) # differential item functioning.

In our meeting on 28/June/2021, we noted that by adding Word Forms, we might be able to model some of the intermediate levels of syntactic proficiency between the 0 and 1 items on the complexity scale. I’m going to start with adding word endings and combines as those seem like the most straight forward. The problem with the word forms and word endings is that we expect them to by non-monotonically related to language proficiency, so I’m no sure how well they’ll fit this model. I’ve also made a new branch on github because I realized I could maintain the ordinal scale of the ordinal items, and use an ordinal IRT model – since binary items are just ordinal items with 2 categories. This appears to be working alright, but I want to keep it in a second branch until I verify this is appropriate (I’m 95% sure it’s right).

Download data from Word Bank:

Inst <- get_instrument_data(language="English (American)", form="WS")
Admin <- get_administration_data(language="English (American)", form="WS")
Item <- get_item_data(language="English (American)", form = "WS")

The American English CDI data has 7 types of items other than words, which are also listed in the code above: combine, complexity, word_endings, word_forms_nouns, word_forms_verbs, word_endings_nouns, and word_endings_verbs. The word_endings type seems to be questions about how often a child applies a particular morphological rule (e.g., adding an s to possessives). The word_forms_verbs and word_forms_nouns seem to be about whether a child can produce specific irregular past tenses or plurals. The word_endings_nouns and word_endings_verbs seem to be about overregularizations on nouns and verbs.

Combine data sets. Create binary variable for accuracy.

Complex <- Admin %>%
  full_join(.,Inst, by="data_id") %>%
  full_join(., Item, by="num_item_id") %>%
  filter(longitudinal==FALSE) %>%
  filter(type == "combine" | 
           type == "complexity" |
           type == "word_endings" |
           type== "word_forms_nouns" |
           type == "word_forms_verbs" | 
           type == "word_endings_nouns" |
           type == "word_endings_verbs"
         ) %>%
  mutate(
    out = ifelse(value=="complex" | value=="sometimes" | value=="produces", yes=1, 
                 no = ifelse(value=="often", yes=2, no =0))
  )

As a quick sanity check, let’s see which items have ordinal data (to make sure the above code worked correctly).

Complex %>%
  filter(value=="often") %>%
  group_by(definition) %>%
  count()

Ok, it’s only the combine variable and word endings variables – that’s right.

0.1 Data Screening

Let’s see how many missing values there are.

Complex %>%
  filter(is.na(out)) %>%
  group_by(definition) %>%
  count()

Looks like there are 1426 participants don’t have item-level complexity scores. There are a further 2 who don’t have word forms, and 1 who doesn’t have word endings.

Complex$complexity_category <- ifelse(Complex$complexity_category == "", yes=Complex$type, no=Complex$complexity_category)

Let’s look at means and SDs for each item:

Complex %>% 
  filter(!is.na(out)) %>%
  group_by(definition) %>%
  summarise(
    mean=mean(out), 
    sd=sd(out),
    sum=sum(out),
    category = first(complexity_category)
  ) %>%
  arrange(category) %>%
  kable(caption="Means and SDs for Each Item (Arranged by Item)")
Means and SDs for Each Item (Arranged by Item)
definition mean sd sum category
has your child begun to combine word yet? 1.2833572 0.7960586 3578 combine
baby blanket / baby’s blanket 0.3873745 0.4872378 1080 morphology
daddy car / daddy’s car 0.4451220 0.4970684 1241 morphology
daddy pick me up / daddy picked me up 0.1743185 0.3794512 486 morphology
doggie kiss me / doggie kissed me 0.1879484 0.3907410 524 morphology
I fall down / I fell down 0.1807747 0.3849005 504 morphology
I make tower / I making tower 0.2360115 0.4247055 658 morphology
kitty go away / kitty went away 0.1280488 0.3342041 357 morphology
kitty sleep / kitty sleeping 0.4192970 0.4935326 1169 morphology
more cookie / more cookies 0.3758967 0.4844404 1048 morphology
these my tooth / these my teeth 0.4533716 0.4979103 1264 morphology
two foot / two feet 0.4214491 0.4938798 1175 morphology
two shoe / two shoes 0.4422525 0.4967431 1233 morphology
baby crying / baby crying cuz she’s sad 0.1599713 0.3666452 446 syntax
baby crying / baby is crying 0.1624821 0.3689586 453 syntax
baby want eat / baby want to eat 0.2650646 0.4414468 739 syntax
coffee hot / that coffee hot 0.1929699 0.3947004 538 syntax
cookie mommy / cookie for mommy 0.2413917 0.4280040 673 syntax
doggie table / doggie on table 0.3511478 0.4774147 979 syntax
don’t read book / don’t want you read that book 0.1818508 0.3857907 507 syntax
go bye-bye / wanna go bye-bye 0.2679340 0.4429625 747 syntax
I like read stories / I like to read stories 0.1779053 0.3825016 496 syntax
I no do it / I can’t do it 0.2482066 0.4320498 692 syntax
I sing song / I sing song for you 0.1449067 0.3520699 404 syntax
I want that / I want that one you got 0.1137016 0.3175054 317 syntax
lookit / lookit what I got 0.2410330 0.4277870 672 syntax
lookit me / lookit me dancing 0.2596844 0.4385400 724 syntax
no wash dolly / don’t wash dolly 0.2116212 0.4085310 590 syntax
read me story Mommy / read me a story Mommy 0.1635581 0.3699405 456 syntax
that my truck / that’s my truck 0.2822812 0.4501902 787 syntax
there a kitty / there’s a kitty 0.3210187 0.4669517 895 syntax
turn on light / turn on light so I can see 0.0961263 0.2948172 268 syntax
want cookies / want cookies and milk 0.2048063 0.4036324 571 syntax
want more juice / want juice in there 0.2051650 0.4038946 572 syntax
we made this / me and Paul made this 0.1603300 0.3669776 447 syntax
where’s my dolly / where’s my dolly name Sam 0.1093974 0.3121932 305 syntax
where mommy go / where did mommy go 0.1822095 0.3860863 508 syntax
you fix it / can you fix it 0.1025825 0.3034672 286 syntax
ed 0.3589470 0.5835291 1009 word_endings
ing 0.6143721 0.7237652 1727 word_endings
splural 0.8701530 0.8247950 2446 word_endings
spossess 0.9558876 0.8877438 2687 word_endings
blockses 0.0211697 0.1439757 59 word_endings_nouns
childrens 0.0312164 0.1739332 87 word_endings_nouns
childs 0.0093290 0.0961526 26 word_endings_nouns
feets 0.1729458 0.3782683 482 word_endings_nouns
foots 0.1094367 0.3122423 305 word_endings_nouns
mans 0.0940079 0.2918921 262 word_endings_nouns
mens 0.0200933 0.1403446 56 word_endings_nouns
mices 0.0114819 0.1065557 32 word_endings_nouns
mouses 0.0606387 0.2387092 169 word_endings_nouns
shoeses 0.0782484 0.2686103 218 word_endings_nouns
sockses 0.0692501 0.2539245 193 word_endings_nouns
teeths 0.1395766 0.3466094 389 word_endings_nouns
toeses 0.0581270 0.2340255 162 word_endings_nouns
tooths 0.0534625 0.2249943 149 word_endings_nouns
ated 0.0200933 0.1403446 56 word_endings_verbs
blewed 0.0118407 0.1081882 33 word_endings_verbs
blowed 0.1115895 0.3149173 311 word_endings_verbs
breaked 0.0696089 0.2545324 194 word_endings_verbs
bringed 0.0487980 0.2154841 136 word_endings_verbs
broked 0.0900610 0.2863205 251 word_endings_verbs
buyed 0.0617151 0.2406806 172 word_endings_verbs
camed 0.0107643 0.1032096 30 word_endings_verbs
comed 0.0186581 0.1353385 52 word_endings_verbs
doed 0.0448511 0.2070141 125 word_endings_verbs
dranked 0.0276283 0.1639347 77 word_endings_verbs
drinked 0.0843201 0.2779171 235 word_endings_verbs
eated 0.0559742 0.2299131 156 word_endings_verbs
falled 0.0961607 0.2948645 268 word_endings_verbs
flied 0.0211697 0.1439757 59 word_endings_verbs
getted 0.0211697 0.1439757 59 word_endings_verbs
goed 0.0487980 0.2154841 136 word_endings_verbs
gotted 0.0125583 0.1113779 35 word_endings_verbs
haved 0.0143524 0.1189598 40 word_endings_verbs
heared 0.0200933 0.1403446 56 word_endings_verbs
holded 0.0362397 0.1869195 101 word_endings_verbs
losed 0.0211697 0.1439757 59 word_endings_verbs
losted 0.0165052 0.1274308 46 word_endings_verbs
maked 0.0423394 0.2013985 118 word_endings_verbs
ranned 0.0193757 0.1378662 54 word_endings_verbs
runned 0.0441335 0.2054285 123 word_endings_verbs
satted 0.0060997 0.0778763 17 word_endings_verbs
seed 0.0279871 0.1649654 78 word_endings_verbs
sitted 0.0258342 0.1586690 72 word_endings_verbs
taked 0.0405454 0.1972699 113 word_endings_verbs
wented 0.0086114 0.0924138 24 word_endings_verbs
children 0.1302943 0.3366874 363 word_forms_nouns
feet 0.5681981 0.4954161 1583 word_forms_nouns
men 0.1162958 0.3206368 324 word_forms_nouns
mice 0.0972721 0.2963811 271 word_forms_nouns
teeth 0.6977746 0.4593048 1944 word_forms_nouns
ate 0.2882268 0.4530185 803 word_forms_verbs
blew 0.1335248 0.3402021 372 word_forms_verbs
bought 0.1166547 0.3210660 325 word_forms_verbs
broke 0.4605169 0.4985281 1283 word_forms_verbs
came 0.1148600 0.3189102 320 word_forms_verbs
drank 0.1765973 0.3813960 492 word_forms_verbs
drove 0.0717875 0.2581821 200 word_forms_verbs
fell 0.3646805 0.4814269 1016 word_forms_verbs
flew 0.1105528 0.3136338 308 word_forms_verbs
got 0.2781766 0.4481813 775 word_forms_verbs
had 0.1256281 0.3314893 350 word_forms_verbs
heard 0.1127064 0.3162903 314 word_forms_verbs
held 0.0678392 0.2515149 189 word_forms_verbs
lost 0.2333094 0.4230134 650 word_forms_verbs
made 0.2132089 0.4096475 594 word_forms_verbs
ran 0.1819813 0.3858984 507 word_forms_verbs
sat 0.1593683 0.3660849 444 word_forms_verbs
saw 0.1870065 0.3899867 521 word_forms_verbs
took 0.1679828 0.3739181 468 word_forms_verbs
went 0.1719311 0.3773883 479 word_forms_verbs

Let’s also check the distribution of means for each category

Complex %>%
  filter(!is.na(out)) %>%
  group_by(definition) %>%
  summarise(
    mean=mean(out), 
    category = first(complexity_category),
    cat2 = str_replace(category, "_", " " ) %>% str_replace("_", "\n" ) 
  ) %>%
  ggplot(aes(x=cat2, y= mean)) + geom_boxplot()

It seems like the word endings and the combine words variable have higher means than the morphological and syntactic complexity items (though they are ordinal, but when I treated these as binary in an older version the results were the same). The over-regularization items have very low means, and the word forms looks a bit like the syntactic and morphological complexity.

One difficulty I just thought of is that the we would expect the word form and over-regularizations to be non-monotonially related to overall grammatical proficiency, given the inverse u shaped curve we typically see in the development of morphology.

Therefore, I’m going to first add the combine item and the word endings variables, see how that changes the results of the IRT model, and then plot the over-regularization and word form scores against this, to see if we see the non-monotonicity.

0.2 Data Preparation

Prepare data set for IRT modeling. I’m going to re-name each item to its category (morphology or syntax) and its item number so as to make some graphs easier to read.

Complex_short_with_ids <- Complex %>% 
  dplyr::select(data_id, value, out, complexity_category, num_item_id) %>%
  mutate(
    label = str_c(complexity_category, num_item_id)
  ) %>%
  pivot_wider(id_cols=data_id, names_from = "label", values_from="out") %>%
  drop_na()


Complex_short <- Complex_short_with_ids %>%
  dplyr::select(starts_with(c("combine", "morphology", "syntax")), c("word_endings686", "word_endings687", "word_endings688", "word_endings689")) # dataset for IRT can't have IDs

0.3 Dimensionality Assessment

Prior to fitting the IRT model, it’s worth looking at the correlation structure of the data to get a sense of the possible dimensionality.

Create tetrachoric correlation matrix (for binary data). I created a vector with labels, for labeling of the correlation plot.

Complex_poly <- polychoric(Complex_short)
## Warning in polychoric(Complex_short): The items do not have an equal number of
## response alternatives, global set to FALSE.
## Warning in matpLower(x, nvar, gminx, gmaxx, gminy, gmaxy): 13 cells were
## adjusted for 0 values using the correction for continuity. Examine your data
## carefully.
## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was done
rho <- Complex_poly$rho

lab = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42") # create labels for corPlot because the defintions are long and take up a ton of space

Check out the correlations

corPlot(rho, labels=lab)

They look really high. Let’s examine the dimensionality, using princals, which is basically just a version of PCA that is appropriate for categorical data.

pc <- princals(rho)

plot(pc)

Seems like there are more intermediate points in between the two clouds of vectors I saw with the earlier version of the data set. However, the plot has fipped, so that most of the vectors are pointing toward the positive first component.

fa.parallel(rho, fa="fa", cor="poly", n.obs = 2785)
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Parallel analysis suggests that the number of factors =  7  and the number of components =  NA

Parallel analysis suggests one very large factor and a small second one. Although with these additional items, I’m seeing some convergence complaints that I wasn’t seeing before. These dimensionality tests aren’t of substantive interest, so I’m not especially worried, but it’s worth keeping this in mind for the later models.

A third metric is very simple structure. This tells us the best solution assuming different values of k, where k is the maximum number of latent variables an observed variable loads on.

vss(rho, cor="poly", n.obs = 2785)
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## 
## Very Simple Structure
## Call: vss(x = rho, n.obs = 2785, cor = "poly")
## VSS complexity 1 achieves a maximimum of 1  with  1  factors
## VSS complexity 2 achieves a maximimum of 1  with  2  factors
## 
## The Velicer MAP achieves a minimum of 0.02  with  3  factors 
## BIC achieves a minimum of  74868.12  with  8  factors
## Sample Size adjusted BIC achieves a minimum of  76625.19  with  8  factors
## 
## Statistics by number of factors 
##   vss1 vss2   map dof  chisq prob sqresid fit RMSEA    BIC  SABIC complex
## 1 1.00 0.00 0.029 819 106986    0    3.38   1  0.22 100489 103092     1.0
## 2 0.50 1.00 0.018 778  96146    0    1.89   1  0.21  89975  92447     1.8
## 3 0.38 0.85 0.016 738  91813    0    1.50   1  0.21  85959  88304     2.4
## 4 0.38 0.83 0.017 699  88809    0    1.30   1  0.21  83265  85486     2.6
## 5 0.36 0.76 0.017 661  85563    0    1.10   1  0.21  80320  82421     3.0
## 6 0.38 0.77 0.017 624  82925    0    0.97   1  0.22  77975  79958     3.0
## 7 0.38 0.77 0.017 588  80796    0    0.85   1  0.22  76132  78000     3.1
## 8 0.38 0.78 0.019 553  79255    0    0.78   1  0.23  74868  76625     3.1
##   eChisq   SRMR eCRMS  eBIC
## 1   4890 0.0319 0.033 -1606
## 2   1903 0.0199 0.021 -4268
## 3   1310 0.0165 0.018 -4544
## 4   1045 0.0148 0.016 -4499
## 5    784 0.0128 0.015 -4459
## 6    631 0.0115 0.013 -4319
## 7    493 0.0101 0.012 -4171
## 8    424 0.0094 0.012 -3963

Very simple structure also suggests one, if we assume each variable loads on one factor only. If we assume everything loads on 2 factors then a 2 factor solution is best. Again, note that the algorithm is having difficulty with estimation.

1 IRT Modeling

1.1 2PL model

Let’s first fit the 2pl model and examine the output.

m1 <- mirt(Complex_short, 1, itemtype="graded", verbose=FALSE) # graded because we are treating all variables as ordinal. 
m1
## 
## Call:
## mirt(data = Complex_short, model = 1, itemtype = "graded", verbose = FALSE)
## 
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 341 EM iterations.
## mirt version: 1.33.2 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 61
## Latent density type: Gaussian 
## 
## Log-likelihood = -34256.68
## Estimated parameters: 89 
## AIC = 68691.36; AICc = 68697.31
## BIC = 69219.31; SABIC = 68936.53
## G2 (1e+10) = 35410.77, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN
summary(m1)
##                    F1    h2
## combine760      0.909 0.826
## morphology761   0.931 0.867
## morphology762   0.899 0.808
## morphology763   0.913 0.833
## morphology764   0.930 0.865
## morphology765   0.940 0.883
## morphology766   0.885 0.783
## morphology767   0.906 0.821
## morphology768   0.912 0.832
## morphology769   0.920 0.846
## morphology770   0.942 0.887
## morphology771   0.955 0.911
## morphology772   0.933 0.870
## syntax773       0.959 0.919
## syntax774       0.938 0.880
## syntax775       0.949 0.901
## syntax776       0.956 0.914
## syntax777       0.954 0.909
## syntax778       0.941 0.886
## syntax779       0.900 0.810
## syntax780       0.952 0.906
## syntax781       0.932 0.869
## syntax782       0.935 0.875
## syntax783       0.957 0.916
## syntax784       0.949 0.900
## syntax785       0.968 0.938
## syntax786       0.968 0.937
## syntax787       0.964 0.930
## syntax788       0.943 0.889
## syntax789       0.938 0.880
## syntax790       0.968 0.937
## syntax791       0.972 0.944
## syntax792       0.947 0.897
## syntax793       0.966 0.932
## syntax794       0.935 0.874
## syntax795       0.915 0.838
## syntax796       0.968 0.937
## syntax797       0.954 0.910
## word_endings686 0.839 0.705
## word_endings687 0.865 0.748
## word_endings688 0.845 0.714
## word_endings689 0.886 0.784
## 
## SS loadings:  36.511 
## Proportion Var:  0.869 
## 
## Factor correlations: 
## 
##    F1
## F1  1

Let’s check the overall fit of the model. The chi square statistic here isn’t super meaningful given the sample size.

M2(m1)

RMSEA and CFI are good, SRMSR is a bit high but pretty reasonable. We can see if individual items misfit – i.e., do the observed cell counts differ from the expected cell counts.

itemfit(m1)

A few items significantly misfit. Let’s take everything with a p value less than .10. We’ll look at item gam plots of those, which plot the item response probability against the latent variable using a gam function, to look for deviations in the item-response function.

misfit <- itemfit(m1) %>% # Get labels of mis-fitting items. 
  filter(p.S_X2 <= .10) %>%
  dplyr::select(item) %>%
  as.vector()

items_good <- dplyr::select(Complex_short, -all_of(misfit$item)) # Well fitting items
items_bad <- dplyr::select(Complex_short, all_of(misfit$item)) # Poorly fitting items

mod_fit <- mirt(items_good, 1, "graded", verbose=FALSE) # Calculate factor scores using only the well fitting items. 
Theta <- fscores(mod_fit)

2 Plot non-linear item response theory curves with the factor scores

IG760 <- itemGAM(items_bad$combine760, Theta)

IG762 <- itemGAM(items_bad$morphology762,Theta)
IG769<- itemGAM(items_bad$morphology769,Theta)
IG771<- itemGAM(items_bad$morphology771,Theta)


IG775 <- itemGAM(items_bad$syntax775,Theta)
IG776 <- itemGAM(items_bad$syntax776,Theta)
IG785 <- itemGAM(items_bad$syntax785,Theta)
IG787<- itemGAM(items_bad$syntax787,Theta)

IG688 <- itemGAM(items_bad$word_endings688,Theta)
IG689<- itemGAM(items_bad$word_endings689,Theta)

Combine

plot(IG760)

Flattens out pretty quickly at a difficulty of 0.

plot(IG762); plot(IG769); plot(IG771)

Syntax items

plot(IG775); plot(IG776); plot(IG785); plot(IG787)

plot(IG688); plot(IG689)

Next let’s look at the relationship between the true scores and the raw scores.

true_raw <- Complex_short %>%
  mutate(
    Raw = rowSums(.[,1:42]), 
    Theta = fscores(m1, method="MAP")
  ) %>%
  ggplot(aes(x=Theta, y=Raw)) + geom_point() + stat_smooth(method="loess") + theme_minimal()

raw_score <- Complex_short %>%
  mutate(
    Raw = rowSums(.[,1:42])
    ) %>% 
   ggplot(aes(x=Raw)) + geom_histogram() + theme_minimal()

true_score <- Complex_short %>%
  mutate(
    Theta = fscores(m1, method="MAP")
    ) %>% 
   ggplot(aes(x=Theta)) + geom_histogram() + theme_minimal()

library(patchwork)

true_raw/(true_score + raw_score) 
## `geom_smooth()` using formula 'y ~ x'
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Well, we’re getting a little more sensitivity at the lower end of the scale than we had before, and the results are starting to look a bit closer to gaussian-ish.

Next let’s take a look at the difficulty and discrimination parameter for each item.

Note to self: This block of code will need to be fixed if I treat variables as ordinal, as the output of coef is now different.

#coefs_2pl <- coef(m1, as.data.frame = TRUE) %>% 
# t() %>%
#  as_tibble() %>%
#  dplyr::select(-c(149:150)) %>%
#  pivot_longer(everything()) %>%
#  separate(, col=name, into=c("item", "parameter"), sep="([.])") %>%
#  pivot_wider(id_cols=item, names_from=parameter, values_from=value) %>%
#  dplyr::select(item, a1, d) %>%
#  mutate(
#   category =  gsub('[[:digit:]]+', '', item)
#  )
  
#ggplot(coefs_2pl,  
#       aes(x = a1, y = d)) + 
#  geom_point(alpha = .3) + 
#  ggrepel::geom_text_repel(data = coefs_2pl, 
#                  aes(label = item), size = 3) + 
#  xlab("Discrimination") + 
#  ylab("Difficulty") + theme_minimal()
#ggplot(coefs_2pl, aes(x=a1, fill=category)) + geom_histogram() + theme_minimal()
#ggplot(coefs_2pl, aes(x=d, fill=category)) + geom_histogram() + theme_minimal()

I was worried that the over-regularization and irregular items would be non-monotonically related to the factor scores, since we expect children to go through qualitatively distinct stages. So I summed these sections and plotted them against the factor score below.

Theta_by_Wordforms <- Complex_short_with_ids %>%
  mutate(
    rn = row_number()
  ) %>%
  transmute(
    Word_form = rowSums(.[,6:30]), 
    Over_Reg = rowSums(.[,31:75]),
    Theta = as.numeric(fscores(m1, "MAP"))
  ) 

p1 <- ggplot(Theta_by_Wordforms, aes(x=Theta, y=Word_form)) + geom_point() + stat_smooth(method="loess")
p2 <- ggplot(Theta_by_Wordforms, aes(x=Theta, y=Over_Reg)) + geom_point() + stat_smooth(method="loess")
p3 <- ggplot(Theta_by_Wordforms, aes(x=Theta)) + geom_histogram()
p4 <- ggplot(Theta_by_Wordforms, aes(x=Word_form)) + geom_histogram()
p5 <- ggplot(Theta_by_Wordforms, aes(x=Over_Reg)) + geom_histogram()
             
(p1 | p2)/(p3 | p4 | p5)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Interesting to see that neither is clearly non-monotonic, even though that’s what we would expect from theories of language acquisition.

Now I’m going to try adding each of those sets of variables to the model I fit above. I do the word forms first and then the over-regularization

3 Adding Word Forms.

Complex_short_with_ids <- Complex %>% 
  dplyr::select(data_id, value, out, complexity_category, num_item_id) %>%
  mutate(
    label = str_c(complexity_category, num_item_id)
  ) %>%
  pivot_wider(id_cols=data_id, names_from = "label", values_from="out") %>%
  drop_na()


Complex_short <- Complex_short_with_ids %>%
  dplyr::select(starts_with(c("combine", "morphology", "syntax", "word_forms")), c("word_endings686", "word_endings687", "word_endings688", "word_endings689")) # dataset for IRT can't have IDs
Complex_poly <- polychoric(Complex_short)
## Warning in polychoric(Complex_short): The items do not have an equal number of
## response alternatives, global set to FALSE.
## Warning in matpLower(x, nvar, gminx, gmaxx, gminy, gmaxy): 13 cells were
## adjusted for 0 values using the correction for continuity. Examine your data
## carefully.
## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was done
rho <- Complex_poly$rho

lab = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60",
"61", "62", "63", "64", "65", "66", "67") # create labels for corPlot because the defintions are long and take up a ton of space
corPlot(rho, labels=lab)

pc <- princals(rho)

plot(pc)

Looks like there’s a big clumb of vectors in one area but then way more pointing out in random directions, which I think will create some fit issues.

fa.parallel(rho, fa="fa", cor="poly", n.obs = 2785)
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Parallel analysis suggests that the number of factors =  8  and the number of components =  NA
vss(rho, fa="fa", cor="tet", n.obs = 2785)
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## 
## Very Simple Structure
## Call: vss(x = rho, n.obs = 2785, cor = "tet", fa = "fa")
## VSS complexity 1 achieves a maximimum of 0.99  with  1  factors
## VSS complexity 2 achieves a maximimum of 1  with  2  factors
## 
## The Velicer MAP achieves a minimum of 0.01  with  5  factors 
## BIC achieves a minimum of  359826.3  with  8  factors
## Sample Size adjusted BIC achieves a minimum of  365237.3  with  8  factors
## 
## Statistics by number of factors 
##   vss1 vss2   map  dof  chisq prob sqresid  fit RMSEA    BIC  SABIC complex
## 1 0.99 0.00 0.035 2144 426325    0    14.4 0.99  0.27 409319 416131     1.0
## 2 0.56 1.00 0.018 2078 403934    0     7.0 1.00  0.26 387451 394054     1.7
## 3 0.37 0.85 0.013 2013 391590    0     4.7 1.00  0.26 375623 382019     2.2
## 4 0.37 0.85 0.013 1949 385741    0     4.1 1.00  0.27 370282 376474     2.3
## 5 0.37 0.84 0.012 1886 382110    0     3.4 1.00  0.27 367150 373142     2.4
## 6 0.38 0.83 0.012 1824 378989    0     3.1 1.00  0.27 364521 370316     2.5
## 7 0.38 0.83 0.012 1763 376197    0     2.8 1.00  0.28 362213 367815     2.5
## 8 0.38 0.83 0.013 1703 373335    0     2.6 1.00  0.28 359826 365237     2.5
##   eChisq  SRMR eCRMS   eBIC
## 1  26667 0.047 0.047   9661
## 2   9604 0.028 0.029  -6879
## 3   5044 0.020 0.021 -10923
## 4   3928 0.018 0.019 -11531
## 5   3050 0.016 0.017 -11909
## 6   2636 0.015 0.016 -11832
## 7   2352 0.014 0.015 -11632
## 8   2065 0.013 0.015 -11443
m2 <- mirt(Complex_short, 1, itemtype="graded", verbose=FALSE)

summary(m2)
##                        F1    h2
## combine760          0.894 0.799
## morphology761       0.924 0.853
## morphology762       0.894 0.800
## morphology763       0.900 0.809
## morphology764       0.919 0.845
## morphology765       0.932 0.869
## morphology766       0.882 0.778
## morphology767       0.897 0.805
## morphology768       0.905 0.820
## morphology769       0.905 0.819
## morphology770       0.941 0.885
## morphology771       0.955 0.911
## morphology772       0.927 0.860
## syntax773           0.952 0.906
## syntax774           0.928 0.861
## syntax775           0.939 0.882
## syntax776           0.953 0.909
## syntax777           0.947 0.898
## syntax778           0.930 0.866
## syntax779           0.891 0.794
## syntax780           0.943 0.889
## syntax781           0.925 0.855
## syntax782           0.929 0.864
## syntax783           0.949 0.901
## syntax784           0.944 0.890
## syntax785           0.963 0.927
## syntax786           0.962 0.926
## syntax787           0.961 0.923
## syntax788           0.933 0.871
## syntax789           0.930 0.864
## syntax790           0.961 0.924
## syntax791           0.964 0.929
## syntax792           0.939 0.881
## syntax793           0.960 0.922
## syntax794           0.928 0.861
## syntax795           0.902 0.814
## syntax796           0.962 0.926
## syntax797           0.950 0.902
## word_forms_nouns690 0.811 0.658
## word_forms_nouns691 0.757 0.573
## word_forms_nouns692 0.736 0.541
## word_forms_nouns693 0.699 0.489
## word_forms_nouns694 0.792 0.627
## word_forms_verbs695 0.878 0.772
## word_forms_verbs696 0.836 0.699
## word_forms_verbs697 0.896 0.803
## word_forms_verbs698 0.821 0.674
## word_forms_verbs699 0.885 0.782
## word_forms_verbs700 0.872 0.761
## word_forms_verbs701 0.866 0.751
## word_forms_verbs702 0.865 0.748
## word_forms_verbs703 0.896 0.803
## word_forms_verbs704 0.880 0.774
## word_forms_verbs705 0.914 0.835
## word_forms_verbs706 0.922 0.850
## word_forms_verbs707 0.905 0.818
## word_forms_verbs708 0.893 0.798
## word_forms_verbs709 0.936 0.877
## word_forms_verbs710 0.830 0.688
## word_forms_verbs711 0.873 0.761
## word_forms_verbs712 0.920 0.846
## word_forms_verbs713 0.877 0.768
## word_forms_verbs714 0.910 0.827
## word_endings686     0.834 0.696
## word_endings687     0.855 0.731
## word_endings688     0.846 0.716
## word_endings689     0.894 0.799
## 
## SS loadings:  54.503 
## Proportion Var:  0.813 
## 
## Factor correlations: 
## 
##    F1
## F1  1
M2(m2)
m2b <- mirt(Complex_short, 2, "graded", verbose=FALSE)
summary(m2b, "oblimin", suppress=.20)
## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##                         F1    F2    h2
## combine760          -0.941    NA 0.823
## morphology761       -0.913    NA 0.862
## morphology762       -0.767    NA 0.802
## morphology763       -0.992    NA 0.832
## morphology764       -0.952    NA 0.859
## morphology765       -0.910    NA 0.876
## morphology766       -0.679 0.218 0.772
## morphology767       -0.868    NA 0.813
## morphology768       -0.856    NA 0.828
## morphology769       -1.044    NA 0.848
## morphology770       -0.733 0.223 0.878
## morphology771       -0.750 0.220 0.905
## morphology772       -0.834    NA 0.858
## syntax773           -0.923    NA 0.914
## syntax774           -0.946    NA 0.873
## syntax775           -0.975    NA 0.896
## syntax776           -0.845    NA 0.907
## syntax777           -0.938    NA 0.904
## syntax778           -1.005    NA 0.884
## syntax779           -0.919    NA 0.804
## syntax780           -0.955    NA 0.901
## syntax781           -0.889    NA 0.862
## syntax782           -0.851    NA 0.864
## syntax783           -0.994    NA 0.914
## syntax784           -0.869    NA 0.893
## syntax785           -0.933    NA 0.933
## syntax786           -0.931    NA 0.932
## syntax787           -0.907    NA 0.926
## syntax788           -0.976    NA 0.885
## syntax789           -0.951    NA 0.876
## syntax790           -0.973    NA 0.935
## syntax791           -0.982    NA 0.941
## syntax792           -0.946    NA 0.893
## syntax793           -0.906    NA 0.927
## syntax794           -0.898    NA 0.865
## syntax795           -0.997    NA 0.835
## syntax796           -0.956    NA 0.934
## syntax797           -0.896    NA 0.906
## word_forms_nouns690 -0.438 0.388 0.644
## word_forms_nouns691     NA 0.714 0.637
## word_forms_nouns692     NA 0.835 0.590
## word_forms_nouns693     NA 0.688 0.513
## word_forms_nouns694     NA 0.682 0.686
## word_forms_verbs695     NA 0.764 0.816
## word_forms_verbs696     NA 0.949 0.770
## word_forms_verbs697     NA 0.733 0.826
## word_forms_verbs698     NA 0.749 0.739
## word_forms_verbs699     NA 0.797 0.814
## word_forms_verbs700     NA 0.921 0.830
## word_forms_verbs701     NA 1.013 0.827
## word_forms_verbs702     NA 0.816 0.816
## word_forms_verbs703     NA 0.801 0.837
## word_forms_verbs704 -0.238 0.676 0.799
## word_forms_verbs705     NA 0.839 0.877
## word_forms_verbs706 -0.207 0.749 0.878
## word_forms_verbs707     NA 0.945 0.882
## word_forms_verbs708     NA 0.774 0.840
## word_forms_verbs709 -0.396 0.571 0.884
## word_forms_verbs710  0.226 1.098 0.817
## word_forms_verbs711     NA 0.947 0.838
## word_forms_verbs712 -0.251 0.703 0.870
## word_forms_verbs713     NA 0.908 0.835
## word_forms_verbs714 -0.260 0.680 0.845
## word_endings686     -0.692    NA 0.695
## word_endings687     -0.793    NA 0.737
## word_endings688     -0.582 0.281 0.708
## word_endings689         NA 0.727 0.778
## 
## Rotated SS loadings:  34.058 17.043 
## 
## Factor correlations: 
## 
##        F1     F2
## F1  1.000 -0.888
## F2 -0.888  1.000

It looks like the word forms are loading consistently on a different factor from the morphosyntactic complexity factor.

Another possibility is that all of those items are guesses–especially since they’re means are all low.

4 Adding Overregularizations.

Complex_short_with_ids <- Complex %>% 
  dplyr::select(data_id, value, out, complexity_category, num_item_id) %>%
  mutate(
    label = str_c(complexity_category, num_item_id)
  ) %>%
  pivot_wider(id_cols=data_id, names_from = "label", values_from="out") %>%
  drop_na()


Complex_short <- Complex_short_with_ids %>%
  dplyr::select(starts_with(c("combine", "morphology", "syntax", "word_endings"))) # dataset for IRT can't have IDs
Complex_poly <- polychoric(Complex_short)
## Warning in polychoric(Complex_short): The items do not have an equal number of
## response alternatives, global set to FALSE.
## Warning in matpLower(x, nvar, gminx, gmaxx, gminy, gmaxy): 31 cells were
## adjusted for 0 values using the correction for continuity. Examine your data
## carefully.
## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was done
rho <- Complex_poly$rho

lab = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60",
"61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80",
"81", "82", "83", "84", "85", "86", "87" ) # create labels for corPlot because the defintions are long and take up a ton of space

corPlot(rho, labels=lab)

pc <- princals(rho)
plot(pc)

fa.parallel(rho, fa="fa", cor="tet", n.obs = 2785)
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Parallel analysis suggests that the number of factors =  15  and the number of components =  NA
vss(rho, fa="fa", cor="tet", n.obs = 2785)
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.

## 
## Very Simple Structure
## Call: vss(x = rho, n.obs = 2785, cor = "tet", fa = "fa")
## VSS complexity 1 achieves a maximimum of 0.93  with  1  factors
## VSS complexity 2 achieves a maximimum of 0.99  with  2  factors
## 
## The Velicer MAP achieves a minimum of 0.01  with  5  factors 
## BIC achieves a minimum of  1230373  with  8  factors
## Sample Size adjusted BIC achieves a minimum of  1240137  with  8  factors
## 
## Statistics by number of factors 
##   vss1 vss2   map  dof   chisq prob sqresid  fit RMSEA     BIC   SABIC complex
## 1 0.93 0.00 0.089 3654 1359800    0     169 0.93  0.37 1330816 1342426     1.0
## 2 0.72 0.99 0.017 3568 1293765    0      26 0.99  0.36 1265464 1276801     1.3
## 3 0.70 0.97 0.015 3483 1281702    0      20 0.99  0.36 1254075 1265142     1.4
## 4 0.69 0.97 0.014 3399 1274046    0      17 0.99  0.37 1247085 1257885     1.5
## 5 0.69 0.96 0.014 3316 1268498    0      15 0.99  0.37 1242196 1252732     1.6
## 6 0.68 0.95 0.014 3234 1263953    0      14 0.99  0.37 1238301 1248576     1.7
## 7 0.69 0.94 0.014 3153 1259157    0      13 0.99  0.38 1234147 1244165     1.7
## 8 0.69 0.94 0.014 3073 1254748    0      11 1.00  0.38 1230373 1240137     1.8
##   eChisq  SRMR eCRMS   eBIC
## 1 407904 0.140 0.142 378921
## 2  39389 0.043 0.045  11087
## 3  27362 0.036 0.038   -265
## 4  21675 0.032 0.034  -5286
## 5  19086 0.030 0.032  -7217
## 6  17125 0.029 0.031  -8527
## 7  15278 0.027 0.029  -9731
## 8  13617 0.026 0.028 -10758
m3b <- mirt(Complex_short, 2, "graded", verbose=FALSE)
## EM cycles terminated after 500 iterations.
summary(m3b, "oblimin", suppress=.20)
## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##                           F1     F2    h2
## combine760                NA  0.978 0.835
## morphology761             NA  0.908 0.867
## morphology762             NA  0.883 0.806
## morphology763             NA  0.954 0.838
## morphology764             NA  0.935 0.865
## morphology765             NA  0.920 0.883
## morphology766             NA  0.901 0.784
## morphology767             NA  0.889 0.822
## morphology768             NA  0.907 0.831
## morphology769             NA  0.940 0.848
## morphology770             NA  0.871 0.893
## morphology771             NA  0.922 0.915
## morphology772             NA  0.914 0.877
## syntax773                 NA  0.965 0.920
## syntax774                 NA  0.943 0.881
## syntax775                 NA  0.971 0.901
## syntax776                 NA  0.946 0.916
## syntax777                 NA  0.984 0.911
## syntax778                 NA  0.974 0.888
## syntax779                 NA  0.854 0.816
## syntax780                 NA  0.944 0.906
## syntax781                 NA  0.910 0.869
## syntax782                 NA  0.950 0.875
## syntax783                 NA  0.936 0.918
## syntax784                 NA  0.920 0.902
## syntax785                 NA  0.965 0.939
## syntax786                 NA  0.919 0.939
## syntax787                 NA  0.897 0.940
## syntax788                 NA  0.871 0.903
## syntax789                 NA  0.907 0.881
## syntax790                 NA  0.944 0.938
## syntax791                 NA  0.955 0.944
## syntax792                 NA  0.916 0.897
## syntax793                 NA  0.937 0.932
## syntax794                 NA  0.901 0.883
## syntax795                 NA  0.880 0.844
## syntax796                 NA  0.928 0.940
## syntax797                 NA  0.901 0.915
## word_endings686           NA  0.782 0.706
## word_endings687           NA  0.865 0.749
## word_endings688           NA  0.718 0.735
## word_endings689       -0.201  0.749 0.820
## word_endings_nouns715 -0.894     NA 0.665
## word_endings_nouns716 -0.602  0.352 0.792
## word_endings_nouns717 -0.776     NA 0.837
## word_endings_nouns718 -0.860     NA 0.664
## word_endings_nouns719 -0.828     NA 0.673
## word_endings_nouns720 -0.552  0.338 0.689
## word_endings_nouns721 -0.551  0.373 0.738
## word_endings_nouns722 -0.660  0.304 0.817
## word_endings_nouns723 -0.538  0.352 0.687
## word_endings_nouns724 -1.046 -0.379 0.667
## word_endings_nouns725 -1.046 -0.356 0.684
## word_endings_nouns726 -0.957 -0.231 0.651
## word_endings_nouns727 -0.986 -0.228 0.700
## word_endings_nouns728 -0.904     NA 0.709
## word_endings_verbs729 -0.903     NA 0.804
## word_endings_verbs730 -0.820     NA 0.786
## word_endings_verbs731 -0.532  0.416 0.776
## word_endings_verbs732 -0.652  0.344 0.866
## word_endings_verbs733 -0.608  0.381 0.850
## word_endings_verbs734 -0.706  0.228 0.784
## word_endings_verbs735 -0.741     NA 0.625
## word_endings_verbs736 -0.834     NA 0.854
## word_endings_verbs737 -0.752     NA 0.787
## word_endings_verbs738 -0.729     NA 0.693
## word_endings_verbs739 -0.782     NA 0.815
## word_endings_verbs740 -0.727  0.217 0.804
## word_endings_verbs741 -0.815     NA 0.840
## word_endings_verbs742 -0.739     NA 0.756
## word_endings_verbs743 -0.546  0.420 0.805
## word_endings_verbs744 -0.899     NA 0.843
## word_endings_verbs745 -0.762     NA 0.730
## word_endings_verbs746 -0.833     NA 0.889
## word_endings_verbs747 -0.758     NA 0.818
## word_endings_verbs748 -0.653  0.341 0.863
## word_endings_verbs749 -0.674  0.288 0.818
## word_endings_verbs750 -0.679  0.325 0.886
## word_endings_verbs751 -0.940     NA 0.843
## word_endings_verbs752 -0.701  0.292 0.873
## word_endings_verbs753 -0.881     NA 0.859
## word_endings_verbs754 -0.633  0.337 0.822
## word_endings_verbs755 -0.778     NA 0.830
## word_endings_verbs756 -0.879     NA 0.916
## word_endings_verbs757 -0.793     NA 0.848
## word_endings_verbs758 -0.662  0.312 0.834
## word_endings_verbs759 -0.833     NA 0.825
## 
## Rotated SS loadings:  27.328 37.616 
## 
## Factor correlations: 
## 
##        F1     F2
## F1  1.000 -0.721
## F2 -0.721  1.000

As was the case with the word form variables, all of the word endings variables loaded onto a separate factor from the morphosyntactic complexity variables.

5 Quick Interim Summary

So a quick re-cap: I was concerned about including the word form and over-regularizations variables to the IRT models because they are theoretically non-monotonically related to syntactic proficiency. We would expect children to initially produce high frequency irregular forms, then go through a stage of over-regularization, and finally to produce the correct irregular word forms again. However, when I plotted the total scores on these subscales against a syntactic proficiency variable derived from the other subscales, I did not see evidence of a non-monotonic relationship. Rather in both cases, subscale scores increased with increasing syntactic proficiency, though it appeared to be leveling off for the over-regularizations.

I then tried fitting two more models, one that included the word form variables and another that included the overgeneralization variables. Exploration of the dimensionality suggested that adding these variables resulted in a much clearer multi-dimensional structure (see categorical PCAs and scree plots). In both cases, exploratory two dimensional IRTs indicated two-factor solutions in which one factor was dominated by the morphosyntactic complexity variables and the other was dominated by the word form or overregularization variables.

It’s difficult to say exactly what this means, as these results could plausibly be due to different psychological demands of learning the relevant structure or some sort of difference in reporting biases across the scales.

6 Final Model

I think I’m going to stick with the model that includes the complexity variable and the word endings variables as they seem to increase the sensitivity of the model to lower proficiency children and do not appear to affect the dimensionality much.

Get the data set back.

Complex_short_with_ids <- Complex %>% 
  dplyr::select(data_id, value, out, complexity_category, num_item_id) %>%
  mutate(
    label = str_c(complexity_category, num_item_id)
  ) %>%
  pivot_wider(id_cols=data_id, names_from = "label", values_from="out") %>%
  drop_na()


Complex_short <- Complex_short_with_ids %>%
  dplyr::select(starts_with(c("combine", "morphology", "syntax")), c("word_endings686", "word_endings687", "word_endings688", "word_endings689")) # dataset for IRT can't have IDs
summary(m1)
##                    F1    h2
## combine760      0.909 0.826
## morphology761   0.931 0.867
## morphology762   0.899 0.808
## morphology763   0.913 0.833
## morphology764   0.930 0.865
## morphology765   0.940 0.883
## morphology766   0.885 0.783
## morphology767   0.906 0.821
## morphology768   0.912 0.832
## morphology769   0.920 0.846
## morphology770   0.942 0.887
## morphology771   0.955 0.911
## morphology772   0.933 0.870
## syntax773       0.959 0.919
## syntax774       0.938 0.880
## syntax775       0.949 0.901
## syntax776       0.956 0.914
## syntax777       0.954 0.909
## syntax778       0.941 0.886
## syntax779       0.900 0.810
## syntax780       0.952 0.906
## syntax781       0.932 0.869
## syntax782       0.935 0.875
## syntax783       0.957 0.916
## syntax784       0.949 0.900
## syntax785       0.968 0.938
## syntax786       0.968 0.937
## syntax787       0.964 0.930
## syntax788       0.943 0.889
## syntax789       0.938 0.880
## syntax790       0.968 0.937
## syntax791       0.972 0.944
## syntax792       0.947 0.897
## syntax793       0.966 0.932
## syntax794       0.935 0.874
## syntax795       0.915 0.838
## syntax796       0.968 0.937
## syntax797       0.954 0.910
## word_endings686 0.839 0.705
## word_endings687 0.865 0.748
## word_endings688 0.845 0.714
## word_endings689 0.886 0.784
## 
## SS loadings:  36.511 
## Proportion Var:  0.869 
## 
## Factor correlations: 
## 
##    F1
## F1  1
M2(m1)

6.1 More formal tests of the dimensionality

m1b <- mirt(Complex_short, 2, "graded", verbose=FALSE)
anova(m1, m1b)
## 
## Model 1: mirt(data = Complex_short, model = 1, itemtype = "graded", verbose = FALSE)
## Model 2: mirt(data = Complex_short, model = 2, itemtype = "graded", verbose = FALSE)
summary(m1b, "oblimin", suppress=.2)
## 
## Rotation:  oblimin 
## 
## Rotated factor loadings: 
## 
##                     F1     F2    h2
## combine760      -0.882     NA 0.822
## morphology761   -0.883     NA 0.876
## morphology762   -0.861     NA 0.802
## morphology763   -0.818  0.389 0.923
## morphology764   -0.894     NA 0.859
## morphology765   -0.920     NA 0.865
## morphology766   -0.858     NA 0.761
## morphology767   -0.861     NA 0.815
## morphology768   -0.894     NA 0.821
## morphology769   -0.843  0.327 0.908
## morphology770   -0.916     NA 0.877
## morphology771   -0.931     NA 0.907
## morphology772   -0.906     NA 0.863
## syntax773       -0.947     NA 0.905
## syntax774       -0.906     NA 0.875
## syntax775       -0.936     NA 0.883
## syntax776       -0.954     NA 0.900
## syntax777       -0.945     NA 0.893
## syntax778       -0.939     NA 0.869
## syntax779       -0.895     NA 0.784
## syntax780       -0.925     NA 0.899
## syntax781       -0.934     NA 0.854
## syntax782       -0.919     NA 0.854
## syntax783       -0.958     NA 0.901
## syntax784       -0.953     NA 0.888
## syntax785       -0.967     NA 0.926
## syntax786       -0.979     NA 0.934
## syntax787       -0.982     NA 0.938
## syntax788       -0.961 -0.211 0.902
## syntax789       -0.942     NA 0.866
## syntax790       -0.973     NA 0.929
## syntax791       -0.972     NA 0.935
## syntax792       -0.951     NA 0.884
## syntax793       -0.973     NA 0.926
## syntax794       -0.939     NA 0.859
## syntax795       -0.922     NA 0.827
## syntax796       -0.977     NA 0.932
## syntax797       -0.962     NA 0.902
## word_endings686 -0.777  0.234 0.717
## word_endings687 -0.785  0.320 0.800
## word_endings688 -0.816     NA 0.690
## word_endings689 -0.863     NA 0.726
## 
## Rotated SS loadings:  35.26 0.852 
## 
## Factor correlations: 
## 
##        F1     F2
## F1  1.000 -0.163
## F2 -0.163  1.000

We have a similar solution to what we had in the first version of this analysis. The factor model fits better but it’s not at all clear what the second factor is adding. Everything loads strongly on the first factor and very little on the second factor.

Comapre with 1 factor graded model with semi-parametric latent variable specification.

m1_nonpar <- mirt(Complex_short, 1, itemtype="graded", dentype = "EH", verbose=FALSE)

M2(m1_nonpar)
Complex_short %>%
  mutate(
    Theta_EH = fscores(m1_nonpar, method="MAP"), 
    Theta_gauss = fscores(m1, method="MAP")
  ) %>%
  ggplot(aes(x=Theta_EH, y=Theta_gauss)) + geom_point() + stat_smooth(method="loess") + theme_minimal()
## `geom_smooth()` using formula 'y ~ x'

m1_nonpar2 <- mirt(Complex_short, 1, itemtype="graded", dentype = "EHW", verbose=FALSE)
## EM cycles terminated after 2000 iterations.
M2(m1_nonpar2)
Complex_short %>%
  mutate(
    Theta_EH = fscores(m1_nonpar2, method="MAP"), 
    Theta_gauss = fscores(m1, method="MAP")
  ) %>%
  ggplot(aes(x=Theta_EH, y=Theta_gauss)) + geom_point() + stat_smooth(method="loess") + theme_minimal()
## Warning: The following factor score estimates failed to converge successfully:
##     647
## `geom_smooth()` using formula 'y ~ x'

These semi-parametric alternatives look pretty similar. I think in both cases, the Gaussian model is a bit easier to interpret and will likely make it easier to work with other functionality.

7 Measurement invariance.

I mentioned in the last meeting that these models assume measusrement invariance – basically the probability of producing an item is the same for two children of the same estimated ability who differ according to some other grouping variable. I have begun to assess this, using differential item functioning technques using some of the predictors in the data set. I should say, I’m out over my skis a bit here at this point.

Create data sets with predictor variables for

Complex_short_with_ids <- Complex %>% 
  dplyr::select(data_id, value, out, complexity_category, num_item_id) %>%
  mutate(
    label = str_c(complexity_category, num_item_id)
  ) %>%
  dplyr::select(data_id, label, out) %>%
  pivot_wider(id_cols=data_id, names_from = "label", values_from="out") %>%
  drop_na()
Complex_short <- Complex_short_with_ids %>%
  dplyr::select(starts_with(c("combine", "morphology", "syntax")), c("word_endings686", "word_endings687", "word_endings688", "word_endings689")) # dataset for IRT can't have IDs
Complex_with_predictors <- Admin %>%
  dplyr::select(data_id, age, ethnicity, sex) %>%
  mutate(
    female = ifelse(sex=="Female", yes=1, no=0),
    age_y = ifelse(age < 25, yes=1, no=0)
  ) %>%
  right_join(., Complex_short_with_ids, by="data_id")
Complex_with_predictors %>%
  summarise(
    age_missing = sum(is.na(age)), 
    ethnicity_missing = sum(is.na(ethnicity)), # Some ethnicity variables missing. 
    sex_missing = sum(is.na(sex))
  )
table(Complex_with_predictors$age)
## 
##  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30 
## 185 133 323 219 172 110 112 136 420 246 142 138 104 118 227
table(Complex_with_predictors$ethnicity)
## 
##    Asian    Black    Other    White Hispanic 
##       67      222       93     2202      131
table(Complex_with_predictors$sex)
## 
## Female   Male  Other 
##   1372   1413      0

Use lordif to look for DIF. I’ll start with sex first. I’ve not included the code here, the lordif function, which I used, does not allow verbose=FALSE so it produces a lot of extra information.

dif
## Call:
## lordif(resp.data = Complex_short2, group = Complex_with_predictors$female, 
##     criterion = "Chisqr", alpha = 0.01)
## 
##   Number of DIF groups: 2 
## 
##   Number of items flagged for DIF: 2 of 42 
## 
##   Items flagged: 32, 41 
## 
##   Number of iterations for purification: 2 of 10 
## 
##   Detection criterion: Chisqr 
## 
##   Threshold: alpha = 0.01 
## 
##    item ncat  chi12  chi13  chi23
## 1     1    3 0.7743 0.6636 0.3904
## 2     2    2 0.9241 0.6157 0.3270
## 3     3    2 0.9512 0.9867 0.8794
## 4     4    2 0.8690 0.9863 0.9846
## 5     5    2 0.8975 0.8291 0.5495
## 6     6    2 0.4033 0.5775 0.5274
## 7     7    2 0.3206 0.5316 0.5985
## 8     8    2 0.6836 0.8488 0.6875
## 9     9    2 0.8184 0.9651 0.8925
## 10   10    2 0.0538 0.0621 0.1752
## 11   11    2 0.8915 0.8355 0.5593
## 12   12    2 0.0629 0.0300 0.0594
## 13   13    2 0.4924 0.7323 0.6968
## 14   14    2 0.1439 0.3214 0.7141
## 15   15    2 0.3602 0.3547 0.2663
## 16   16    2 0.3432 0.1444 0.0847
## 17   17    2 0.8740 0.1595 0.0562
## 18   18    2 0.4376 0.6299 0.5705
## 19   19    2 0.0655 0.0956 0.2536
## 20   20    2 0.4993 0.7945 0.9514
## 21   21    2 0.9250 0.3804 0.1654
## 22   22    2 0.7012 0.7687 0.5381
## 23   23    2 0.1270 0.2772 0.6261
## 24   24    2 0.3017 0.3705 0.3377
## 25   25    2 0.6176 0.8807 0.9455
## 26   26    2 0.1887 0.3912 0.6995
## 27   27    2 0.8309 0.5944 0.3186
## 28   28    2 0.9316 0.9880 0.8970
## 29   29    2 0.8823 0.9135 0.6901
## 30   30    2 0.9903 0.4593 0.2122
## 31   31    2 0.5563 0.7444 0.6211
## 32   32    2 0.2245 0.0001 0.0000
## 33   33    2 0.1335 0.2428 0.4463
## 34   34    2 0.5069 0.5350 0.3679
## 35   35    2 0.0138 0.0482 0.9686
## 36   36    2 0.7194 0.8807 0.7236
## 37   37    2 0.8068 0.9214 0.7473
## 38   38    2 0.1434 0.2160 0.3364
## 39   39    3 0.8380 0.9265 0.7391
## 40   40    3 0.1360 0.2886 0.6084
## 41   41    3 0.0046 0.0055 0.1243
## 42   42    3 0.9052 0.9743 0.8456
plot(dif)

Only two items have DIF and they look pretty similar across the two groups te be honest.

I’ve turned off the code I’ll median-split the sample by age, and look for DIF there

Every single item has Dif:

dif
## Call:
## lordif(resp.data = Complex_short2, group = Complex_with_predictors$age_y, 
##     criterion = "Chisqr", alpha = 0.01)
## 
##   Number of DIF groups: 2 
## 
##   Number of items flagged for DIF: 42 of 42 
## 
##   Items flagged: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42 
## 
##   Number of iterations for purification: 9 of 10 
## 
##   Detection criterion: Chisqr 
## 
##   Threshold: alpha = 0.01 
## 
##    item ncat  chi12  chi13  chi23
## 1     1    3 0.0000 0.0000 0.0000
## 2     2    2 0.0000 0.0000 0.0000
## 3     3    2 0.0000 0.0000 0.0000
## 4     4    2 0.0000 0.0000 0.0000
## 5     5    2 0.0000 0.0000 0.0000
## 6     6    2 0.0000 0.0000 0.0000
## 7     7    2 0.0000 0.0000 0.0000
## 8     8    2 0.0000 0.0000 0.0000
## 9     9    2 0.0000 0.0000 0.0000
## 10   10    2 0.0000 0.0000 0.0000
## 11   11    2 0.0000 0.0000 0.0000
## 12   12    2 0.0000 0.0000 0.0000
## 13   13    2 0.0091 0.0007 0.0055
## 14   14    2 0.0000 0.0000 0.0000
## 15   15    2 0.0000 0.0000 0.0000
## 16   16    2 0.0009 0.0009 0.0791
## 17   17    2 0.0043 0.0160 0.7610
## 18   18    2 0.0000 0.0000 0.0000
## 19   19    2 0.0000 0.0000 0.0000
## 20   20    2 0.0000 0.0000 0.0000
## 21   21    2 0.0000 0.0000 0.0000
## 22   22    2 0.0000 0.0000 0.0000
## 23   23    2 0.0000 0.0000 0.0008
## 24   24    2 0.0000 0.0000 0.0002
## 25   25    2 0.0000 0.0000 0.0000
## 26   26    2 0.0000 0.0000 0.4918
## 27   27    2 0.0000 0.0000 0.0001
## 28   28    2 0.0000 0.0000 0.0242
## 29   29    2 0.0000 0.0000 0.0006
## 30   30    2 0.0000 0.0000 0.0000
## 31   31    2 0.0000 0.0000 0.0003
## 32   32    2 0.0000 0.0000 0.0000
## 33   33    2 0.0000 0.0000 0.0000
## 34   34    2 0.0000 0.0000 0.0000
## 35   35    2 0.0064 0.0001 0.0011
## 36   36    2 0.0000 0.0000 0.0002
## 37   37    2 0.0000 0.0000 0.0000
## 38   38    2 0.0000 0.0000 0.0001
## 39   39    3 0.0000 0.0000 0.0000
## 40   40    3 0.0000 0.0000 0.0000
## 41   41    3 0.0000 0.0000 0.0000
## 42   42    3 0.0000 0.0000 0.0000

I can’t get plots I had above, because every item has dif – the factor scores can’t be equated in any reasonable way.

Plot the distribution of scores by age group to see what’s going on.

Complex_with_predictors %>%
  dplyr::select(starts_with(c("morphology", "syntax", "combine")) , c("word_endings686", "word_endings687", "word_endings688", "word_endings689",  "age_y")) %>%
  mutate(
    mean = rowSums(.[1:42])
  ) %>%
  ggplot(aes(y=mean, x=as.factor(age_y))) + geom_boxplot()

There are a lot of 0s in the low vocabulary group. I suspect this makes the IRT model difficult to estimate for them, leading to the DIF issues aboe.

8 Plot Theta against age and vocabulary.

by_Participant <-Complex_short_with_ids %>% 
  mutate(
    Theta = fscores(m1, method="MAP")[,1],
    Raw = rowSums(.[2:38])
    ) %>%
  dplyr::select(data_id, Raw, Theta) %>%
  full_join(Complex) %>%
  group_by(data_id) %>%
  slice(1)
## Joining, by = "data_id"
theta_by_age <- by_Participant %>%
  ggplot(aes(x=age, y=Theta)) + geom_point() + stat_smooth(method="loess")

raw_by_age <- by_Participant %>%
  ggplot(aes(x=age, y=Raw)) + geom_point() + stat_smooth(method="loess")

theta_by_age | raw_by_age
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1429 rows containing non-finite values (stat_smooth).
## Warning: Removed 1429 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1429 rows containing non-finite values (stat_smooth).

## Warning: Removed 1429 rows containing missing values (geom_point).

theta_by_prod <- by_Participant %>%
  ggplot(aes(x=production, y=Theta)) + geom_point() + stat_smooth(method="loess")

raw_by_prod <- by_Participant %>%
  ggplot(aes(x=production, y=Raw)) + geom_point() + stat_smooth(method="loess")

theta_by_prod | raw_by_prod
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1429 rows containing non-finite values (stat_smooth).
## Warning: Removed 1429 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1429 rows containing non-finite values (stat_smooth).

## Warning: Removed 1429 rows containing missing values (geom_point).

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Dutch_Netherlands.1252  LC_CTYPE=Dutch_Netherlands.1252   
## [3] LC_MONETARY=Dutch_Netherlands.1252 LC_NUMERIC=C                      
## [5] LC_TIME=Dutch_Netherlands.1252    
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] lordif_0.3-3    rms_6.2-0       SparseM_1.81    Hmisc_4.5-0    
##  [5] Formula_1.2-4   survival_3.1-12 GGally_2.1.1    patchwork_1.1.1
##  [9] knitr_1.33      Gifi_0.3-9      psych_2.1.3     ltm_1.1-1      
## [13] polycor_0.7-10  msm_1.6.8       MASS_7.3-51.6   mirt_1.33.2    
## [17] lattice_0.20-41 forcats_0.5.1   stringr_1.4.0   dplyr_1.0.5    
## [21] purrr_0.3.4     readr_1.4.0     tidyr_1.1.3     tibble_3.1.1   
## [25] ggplot2_3.3.3   tidyverse_1.3.1 wordbankr_0.3.1
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.0-10        colorspace_2.0-0      ellipsis_0.3.1       
##  [4] htmlTable_2.2.1       base64enc_0.1-3       fs_1.5.0             
##  [7] rstudioapi_0.13       farver_2.1.0          Deriv_4.1.3          
## [10] MatrixModels_0.5-0    fansi_0.4.2           mvtnorm_1.1-1        
## [13] lubridate_1.7.10      xml2_1.3.2            codetools_0.2-16     
## [16] splines_4.0.2         mnormt_2.0.2          jsonlite_1.7.2       
## [19] broom_0.7.6           cluster_2.1.0         dbplyr_2.1.1         
## [22] png_0.1-7             compiler_4.0.2        httr_1.4.2           
## [25] backports_1.2.1       assertthat_0.2.1      Matrix_1.2-18        
## [28] cli_2.5.0             htmltools_0.5.1.1     quantreg_5.85        
## [31] tools_4.0.2           gtable_0.3.0          glue_1.4.2           
## [34] Rcpp_1.0.6            cellranger_1.1.0      jquerylib_0.1.4      
## [37] vctrs_0.3.7           nlme_3.1-148          conquer_1.0.2        
## [40] xfun_0.22             rvest_1.0.0           lifecycle_1.0.0      
## [43] dcurver_0.9.2         polspline_1.1.19      zoo_1.8-9            
## [46] scales_1.1.1          hms_1.0.0             sandwich_3.0-0       
## [49] parallel_4.0.2        expm_0.999-6          RColorBrewer_1.1-2   
## [52] yaml_2.2.1            gridExtra_2.3         sass_0.3.1           
## [55] rpart_4.1-15          reshape_0.8.8         latticeExtra_0.6-29  
## [58] stringi_1.5.3         highr_0.9             RMySQL_0.10.21       
## [61] checkmate_2.0.0       permute_0.9-5         rlang_0.4.11         
## [64] pkgconfig_2.0.3       matrixStats_0.58.0    evaluate_0.14        
## [67] labeling_0.4.2        htmlwidgets_1.5.3     quantregGrowth_1.2-1 
## [70] tidyselect_1.1.1      plyr_1.8.6            magrittr_2.0.1       
## [73] R6_2.5.0              generics_0.1.0        multcomp_1.4-17      
## [76] DBI_1.1.1             pillar_1.6.0          haven_2.4.1          
## [79] foreign_0.8-81        withr_2.4.2           mgcv_1.8-36          
## [82] nnet_7.3-14           modelr_0.1.8          crayon_1.4.1         
## [85] utf8_1.2.1            tmvnsim_1.0-2         rmarkdown_2.7        
## [88] jpeg_0.1-8.1          grid_4.0.2            readxl_1.3.1         
## [91] data.table_1.14.0     vegan_2.5-7           reprex_2.0.0         
## [94] digest_0.6.27         GPArotation_2014.11-1 munsell_0.5.0        
## [97] bslib_0.2.4