Hi,
I am trying to make manual code for the Fixed/Random Effects model estimation to understand how the process works.
However, my manual codes for both models never reject the null for the fake data set.
I want to ask for people’s advice here.
My codes are as follows (the first one is the Random effects model, and then the fixed effects model).
SubjectID, DecisionID, Constant, Treatment, SessionID, SwitchPoint, Gender, Complexity = environment();
n_i=length(unique(SubjectID));n_t=length(unique(DecisionID));n_s=length(unique(SessionID));n_i_per_s=Int64(n_i*n_t/n_s)
ObservationID = repeat([1:1:90;],n_s); # needed to fit into the panel data setting
X = [Constant Treatment SessionID SwitchPoint Gender Complexity];
β = [1.0; 0.75; 0.0; 0.01; 0.0; -0.1]; n_β = length(β);
β̂_Mat = zeros(n_β,rep_time);
p_val_Mat_RE = zeros(n_β,rep_time);
for s = 1:rep_time
ϵ = rand(Normal(0,7.5),n_i*n_t);u = repeat(rand(Normal(0,7.5),n_s),inner=n_i_per_s);
Choice = X * β + ϵ + u; df=(;Choice,Treatment,SessionID, SwitchPoint, Gender, Complexity,SubjectID,ObservationID);
lm_POLS = reg(df, @formula(Choice~Treatment+SessionID+SwitchPoint+Gender+Complexity));
v̂ = residuals(lm_POLS,df);
σ̂_v = sqrt(sum(v̂.^2)/(n_i*n_t-n_β));
sum_of_cross_product = 0;
for i = 1:n_s, j = 1:(n_i_per_s-1), k = (j+1):n_i_per_s
sum_of_cross_product = sum_of_cross_product+v̂[n_i_per_s*(i-1)+j]*v̂[n_i_per_s*(i-1)+k]
end
σ̂_u = sqrt(sum_of_cross_product/(n_s*n_i_per_s*(n_i_per_s-1)/2-n_β))
σ̂_ϵ = sqrt(σ̂_v^2 - σ̂_u^2)
Ω̂_s = σ̂_u^2 .* ones(n_i_per_s,n_i_per_s) .+ Diagonal(fill(σ̂_ϵ^2,n_i_per_s)); # covariance structure for each session
Ω̂ = kron(Diagonal(ones(n_s)),Ω̂_s); # covariance structure for the whole session
β̂_RE = inv(X' * inv(Ω̂ ) * X) * X' * inv(Ω̂ ) * Choice;
resid_RE = Choice .- X * β̂_RE;
V̂_RE = inv(X' * inv(Ω̂ ) * X) .* σ̂_ϵ^2;
t_val_vec_RE = β̂_RE ./ sqrt.(diag(V̂_RE))
for l = 1:n_β
p_val_Mat_RE[l,s] = ifelse(t_val_vec_RE[l]>0, 1-pnorm_std(t_val_vec_RE[l])+pnorm_std(-t_val_vec_RE[l]), pnorm_std(t_val_vec_RE[l])+1-pnorm_std(-t_val_vec_RE[l]));
end
β̂_Mat[:,s] = β̂_RE
end
bias = median(β̂_Mat, dims = 2) .- β;
power_β̂_RE = sum(p_val_Mat_RE.<0.05, dims=2)./rep_time
return β̂_Mat, bias, power_β̂_RE
end
@time β̂_Mat3c, bias3c, power3c =Third_c_manual(100)
println(DataFrame(bias = vec(bias3c), power = vec(power3c)))
SubjectID, DecisionID, Constant, Treatment, SessionID, SwitchPoint, Gender, Complexity = environment();
n_i=length(unique(SubjectID));n_t=length(unique(DecisionID));n_s=length(unique(SessionID));n_i_per_s=Int64(n_i*n_t/n_s)
ObservationID = repeat([1:1:90;],n_s); # needed to fit into the panel data setting
X = [Constant Treatment SessionID SwitchPoint Gender Complexity];
β = [1.0; 0.75; 0.0; 0.01; 0.0; -0.1]; n_β = length(β);
β̂_Mat = zeros(n_β,rep_time);
p_val_Mat_FE = zeros(n_β,rep_time);
for s = 1:rep_time
ϵ = rand(Normal(0,7.5),n_i*n_t);u = repeat(rand(Normal(0,7.5),n_s),inner=n_i_per_s);
Choice = X * β + ϵ + u;
I_s = ones(n_i_per_s)
M_s = Diagonal(I_s) .- I_s * I_s' ./ (I_s'*I_s) # demean operator for a session
M = kron(Diagonal(ones(n_s)),M_s) # demean operator for all sessions
A = Symmetric(X' * M * X)
β̂_FE = inv(A) * X' * M * Choice
resid_FE = Choice .- X * β̂_FE
part_resid_FE = collect(Iterators.partition(resid_FE,n_i_per_s)); # partitioning the residuals wrt clusters
part_X = [X[i:min(i+n_i_per_s-1,n_i*n_t),:] for i in 1:n_i_per_s:(n_i*n_t)]; # partitioning the covariates wrt clusters
Ω̂ = zeros(n_β,n_β);
for j = 1:n_s
X_g = part_X[j]
resid_g = part_resid_FE[j]
Ω̂ = Ω̂ + X_g' * resid_g * resid_g' * X_g
end
V̂_FE = inv(A) * Ω̂ * inv(A)
β̂_Mat[:,s] = β̂_FE
t_val_vec_FE = β̂_FE ./ sqrt.(diag(V̂_FE))
for l = 1:n_β
p_val_Mat_FE[l,s] = ifelse(t_val_vec_FE[l]>0, 1-pnorm_std(t_val_vec_FE[l])+pnorm_std(-t_val_vec_FE[l]), pnorm_std(t_val_vec_FE[l])+1-pnorm_std(-t_val_vec_FE[l]));
end
end
bias = median(β̂_Mat, dims = 2) .- β;
power_β̂_FE = sum(p_val_Mat_FE.<0.05, dims=2)./rep_time
return β̂_Mat, bias, power_β̂_FE
end
@time β̂_Mat3d_FEpack, bias3d_FEpack, power3d_FEpack =Third_d_manual(100)
println(DataFrame(bias = vec(bias3d_FEpack), power = vec(power3d_FEpack)))