library(pacman)
p_load("tidyverse",'nlWaldTest','sandwich', 'Rcpp')
popsize<-100
treatments<-c(0, 1)
treatmentvc<-rep(treatments, each=popsize)
ids<-seq(1:(length(treatments)*popsize))
bins<-c(1,0)
q1a<-sample(bins, popsize,replace=T, prob=c(0.4,0.6))
q2a<-sample(bins, popsize,replace=T, prob=c(0.4,0.6))
q3a<-sample(bins, popsize,replace=T, prob=c(0.4,0.6))
q1b<-sample(bins, popsize,replace=T, prob=c(0.3,0.7))
q2b<-sample(bins, popsize,replace=T, prob=c(0.7,0.3))
q3b<-sample(bins, popsize,replace=T, prob=c(0.4,0.5))
q1<-c(q1a,q1b)
q2<-c(q2a,q2b)
q3<-c(q3a,q3b)
df<-tibble::as_tibble(data.frame(ids, treatmentvc, q1,q2,q3))
long_df<-df %>% pivot_longer(cols=c(q1,q2,q3), values_to='Y', names_to = 'Q' ) %>%
rename(T=treatmentvc, id=ids)
long_df
## # A tibble: 600 x 4
## id T Q Y
## <int> <dbl> <chr> <dbl>
## 1 1 0 q1 1
## 2 1 0 q2 1
## 3 1 0 q3 0
## 4 2 0 q1 0
## 5 2 0 q2 1
## 6 2 0 q3 0
## 7 3 0 q1 0
## 8 3 0 q2 0
## 9 3 0 q3 0
## 10 4 0 q1 0
## # … with 590 more rows
fit <- lm(Y ~ 0 + Q + Q:T, data = long_df)
summary(fit)
##
## Call:
## lm(formula = Y ~ 0 + Q + Q:T, data = long_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74 -0.37 -0.30 0.59 0.70
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Qq1 3.700e-01 4.756e-02 7.780 3.24e-14 ***
## Qq2 4.100e-01 4.756e-02 8.621 < 2e-16 ***
## Qq3 3.700e-01 4.756e-02 7.780 3.24e-14 ***
## Qq1:T -7.000e-02 6.726e-02 -1.041 0.298
## Qq2:T 3.300e-01 6.726e-02 4.906 1.20e-06 ***
## Qq3:T 1.178e-16 6.726e-02 0.000 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4756 on 594 degrees of freedom
## Multiple R-squared: 0.4752, Adjusted R-squared: 0.4699
## F-statistic: 89.63 on 6 and 594 DF, p-value: < 2.2e-16
nlWaldTest::nlWaldtest(fit, Vcov = sandwich::vcovHC(fit), df2 = TRUE,
texts = "(b[1]-(b[1]+b[2]+b[3])/3)^2+
(b[2]-(b[1]+b[2]+b[3])/3)^2+
(b[3]-(b[1]+b[2]+b[3])/3)^2;
(b[1]+b[4]-(b[1]+b[2]+b[3]+b[4]+b[5]+b[6])/3)^2+
(b[2]+b[5]-(b[1]+b[2]+b[3]+b[4]+b[5]+b[6])/3)^2+
(b[3]+b[6]-(b[1]+b[2]+b[3]+b[4]+b[5]+b[6])/3)^2")
##
## Wald F test of restrictions on model parameters
##
## data: fit
## F = 6.883, df1 = 2, df2 = 594, p-value = 0.001109