# My manual Fixed/Random Effects Model perform weirdly

**URL:** <https://discourse.julialang.org/t/my-manual-fixed-random-effects-model-perform-weirdly/107705>\
**Category:** New to Julia\
**Created:** [December 16, 2023, 10:13am UTC](https://discourse.julialang.org/t/my-manual-fixed-random-effects-model-perform-weirdly/107705 "2023-12-16T10:13:53Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Wooheon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wooheon/32/206569_2.png) [@Wooheon](https://discourse.julialang.org/u/Wooheon)\
**Post date:** [December 16, 2023, 10:13am UTC](https://discourse.julialang.org/t/my-manual-fixed-random-effects-model-perform-weirdly/107705/1 "2023-12-16T10:13:53Z")

</div>

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).

```function
    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)))

```

```function
    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)))

```