This is a problem of collinearity among your variables.
The lm command automatically places NAs in the beta vector for variables that were not estimated due to colinearity, but PLM does not.
When you have LoadArea*DischargeArea PLM will three variables to your model:
LoadArea + DischargeArea + LoadArea:DischargeArea
After that PLM will demean them.
In this case, and without further information on your data my guess is that one of these variables is perfectly collinear with one of the factors levels in:
as.factor(Laycan.Day.Diff)
In your case I would try to estimate the model without the factor. If it works you know the factors are causing the problem. If it comes to that you can then convert each factor to a explicit 0/1 dummy and add them one by one until you understand where the problem is coming from.
To determine which variables are collinear you could try something like:
require(data.table)
tmp <- data.table(var1=1:10,var2=55:64,userid=rep(c(1,2),5))
cols <- c('var1','var2')
newnames <- c('demeaned_var1','demeaned_var2')
tmp[,(newnames):=.SD-lapply(.SD,mean),.SDcols=cols,by=userid]
cor(tmp[,newnames,with=F])
Line 5 is the demeaning. This other stack overflow post describes the operations of the data.table that i used above in detail.
The output of the code above will be:
>
demeaned_var1 demeaned_var2
demeaned_var1 1 1
demeaned_var2 1 1
This will tell you which demeaned vars are perfectly collinear.