TL/DR: I'm trying to convert a model from SPSS to R and I'm getting different results (Estimates of Fixed Effects), and I haven't been able to figure out what I'm doing wrong. Syntax, code, and results at the end.
For my thesis, I am studying the effect of salt concentration in cooking water on pasta quality using three different types of pastas. So I have two factors, each with three levels (SALT concentration: 0%, 1%, 3%; Pasta TYPE: CP, RL, WT).
Only one pasta type and one salt concentration could be tested at a time. Each pot of pasta is therefore a BATCH. Two batches were prepared for each combination of factor levels (2 * 9).
A random selection of pieces from each batch were tested for quality. These are all recorded separately.
I assume the samples WITHIN each batch will be correlated and would like to account for this. The variation BETWEEN batches is random (i.e. the first batch for WT-1% should not be correlated with the first batch for CP-3%). Because of this, I have coded each batch individually (1-18) instead of 1-2.
I figured out the model using SPSS, but I would like to switch over to R because I have 10 response variables to do this with, and I also have to deal with back-transforming data.
Data Structure
> str(tpa_2)
'data.frame': 115 obs. of 6 variables:
$ Type : Factor w/ 4 levels "CP","CP2","RL",..: 1 1 1 2 2 2 2 2 2 2 ...
$ Batch : Factor w/ 18 levels "1","2","3","4",..: 1 1 1 2 2 2 2 2 2 2 ...
$ Salt : Factor w/ 3 levels "0","1","3": 1 1 1 1 1 1 1 1 1 1 ...
$ Replicate: Factor w/ 2 levels "1","2": 1 1 1 2 2 2 2 2 2 2 ...
$ bn2S : num -1.341 -0.502 -1.448 -1.02 -1.335 ...
$ SampleID : Factor w/ 115 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
SPSS Model
MIXED bn2S BY SaltCode TypeCode
/CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(1000) MXSTEP(10) SCORING(1)
SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0.00001, ABSOLUTE) PCONVERGE(0, ABSOLUTE)
/FIXED = SaltCode TypeCode SaltCode*TypeCode | SSTYPE(3)
/METHOD = REML
/PRINT = CORB SOLUTION
/RANDOM = INTERCEPT | SUBJECT(Batch) COVTYPE(VC).
R Model
library(afex)
m2<-mixed(bn2S ~Salt*Type + (1|Batch), data = tpa_2)
sum.m2<-summary(m2)
sum.m2
- I get the same results with afex::mixed, lme4::lmer, and nlme::lme. The Anova Table wrapper from afex is structured like the Fixed Effects Tests table from SPSS.
Now onto the results.
I get significantly different results between SPSS and R.
SPSS Results
Type III Tests of Fixed Effects
| Source | Numerator df | Denominator df | F | Sig. |
|---|---|---|---|---|
| Intercept | 1 | 3.594 | 18.178 | .016 |
| SaltCode | 2 | 3.310 | 3.974 | .132 |
| TypeCode | 3 | 2.641 | 272.032 | <.001 |
| SaltCode * TypeCode | 6 | 2.348 | 4.217 | .174 |
Estimates of Fixed Effects a
| Parameter | Estimate | Std. Error | df | t | Sig. | 95% CI Lower Bound | 95% CI Upper Bound |
|---|---|---|---|---|---|---|---|
| Intercept | .877 | .100 | 3.616 | 8.790 | .001 | .588 | 1.165 |
| [SaltCode=1] | .591 | .156 | 5.187 | 3.785 | .012 | .194 | .988 |
| [SaltCode=2] | .548 | .151 | 4.750 | 3.628 | .017 | .154 | .943 |
| [SaltCode=3] | 0b | 0 | . | . | . | . | . |
| [TypeCode=1] | -1.965 | .165 | 3.017 | -11.891 | .001 | -2.489 | -1.441 |
| [TypeCode=2] | -1.827 | .151 | 2.077 | -12.114 | .006 | -2.454 | -1.201 |
| [TypeCode=3] | -.817 | .130 | 2.604 | -6.279 | .012 | -1.270 | -.365 |
| [TypeCode=4] | 0b | 0 | . | . | . | . | . |
| [SaltCode=1] * [TypeCode=1] | -.599 | .282 | 6.340 | -2.128 | .075 | -1.280 | .081 |
| [SaltCode=1] * [TypeCode=2] | -.771 | .226 | 2.593 | -3.410 | .052 | -1.560 | .017 |
| [SaltCode=1] * [TypeCode=3] | -.370 | .208 | 3.902 | -1.779 | .152 | -.953 | .213 |
| [SaltCode=1] * [TypeCode=4] | 0b | 0 | . | . | . | . | . |
| [SaltCode=2] * [TypeCode=1] | -.590 | .245 | 3.657 | -2.407 | .080 | -1.297 | .116 |
| [SaltCode=2] * [TypeCode=2] | -.747 | .223 | 2.474 | -3.352 | .058 | -1.550 | .056 |
| [SaltCode=2] * [TypeCode=3] | .019 | .196 | 3.369 | .097 | .928 | -.569 | .607 |
| [SaltCode=2] * [TypeCode=4] | 0b | 0 | . | . | . | . | . |
| [SaltCode=3] * [TypeCode=1] | 0b | 0 | . | . | . | . | . |
| [SaltCode=3] * [TypeCode=2] | 0b | 0 | . | . | . | . | . |
| [SaltCode=3] * [TypeCode=3] | 0b | 0 | . | . | . | . | . |
| [SaltCode=3] * [TypeCode=4] | 0b | 0 | . | . | . | . | . |
b. This parameter is set to zero because it is redundant.
R Results
nice(m2)
Mixed Model Anova Table (Type 3 tests, S-method) Model: bn2S ~ Salt * Type + (1 | Batch) Data: tpa_2
| Effect | df | F | p.value |
|---|---|---|---|
| Salt | 2, 3.02 | 3.97 | .143 |
| Type | 3, 2.93 | 272.03 | <.001 |
| Salt:Type | 6, 2.70 | 4.22 | .150 |
summary(m2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest] Formula: bn2S ~ Salt * Type + (1 | Batch) Data: data
REML criterion at convergence: 105.3
Scaled residuals:
| Min | 1Q | Median | 3Q | Max |
|---|---|---|---|---|
| -3.2729 | -0.6207 | 0.0923 | 0.6020 | 2.5393 |
Random effects:
| Groups | Name | Variance | Std.Dev. |
|---|---|---|---|
| Batch | (Intercept) | 0.002186 | 0.04675 |
| Residual | 0.106238 | 0.32594 |
Number of obs: 115, groups: Batch, 18
Fixed effects:
| Factor | Estimate | Std. Error | df | t value | Pr(>t) |
|---|---|---|---|---|---|
| (Intercept) | -0.15094 | 0.03540 | 3.59414 | -4.264 | 0.01633 |
| Salt1 | 0.03086 | 0.05363 | 4.66163 | 0.575 | 0.59167 |
| Salt2 | 0.09404 | 0.04908 | 3.34529 | 1.916 | 0.14169 |
| Type1 | -0.95426 | 0.07357 | 5.09880 | -12.972 | 4.24e-05 |
| Type2 | -0.92612 | 0.05933 | 2.11691 | -15.611 | 0.00320 |
| Type3 | 0.47305 | 0.05277 | 3.12594 | 8.965 | 0.00249 |
| Salt1:Type1 | -0.02260 | 0.11533 | 7.69583 | -0.196 | 0.84970 |
| Salt2:Type1 | -0.11903 | 0.09941 | 4.24753 | -1.197 | 0.29371 |
| Salt1:Type2 | -0.08504 | 0.08645 | 2.37978 | -0.984 | 0.41446 |
| Salt2:Type2 | -0.16617 | 0.08371 | 2.10267 | -1.985 | 0.17928 |
| Salt1:Type3 | -0.07265 | 0.07992 | 3.89972 | -0.909 | 0.41597 |
| Salt2:Type3 | 0.21058 | 0.07339 | 3.00943 | 2.869 | 0.06385 |
