title: Power Analysis: Comparing GLM to binomial proportion test

text: I am trying to run a power analysis using simulated data to compare the power of a glm versus a binomial proportion test to detect differences in proportions. For example, suppose you have some proportion that decreases by some amount over X number of time steps. .4,.39,.38,.37 . . . . .01 and simulate data over those time steps based on the decrease. Does a glm approach (differences in slope) “outperform” an approach whereby you simply look at proportional differences. I ran simulations and basically summed up the number of times p<.05 divided by the number of trials. I was expecting that a glm approach would have more power because it would be utilizing all the data across multiple time steps, whereas a binomial proportion test is only comparing two populations ( the beginning and the end points). However the results indicated that a binomial proportion test had more power relative to a glm at fewer time steps and that, depending on the simulated decrease in proportions, the relationship between glm power and binomial proportion test power changed. Interestingly, at greater decreases in the proportion, a binomial proportion test seems have greater power compared to a glm, which to me is counterintuitive since the slope of the glm should be greater. I am attaching the code below my questions in case anyone is interested. My questions are: (1) Am I interpreting the results and p-values correctly? (2) If I am interested in trends, does glm results really have lower power and, if so, is there a way to combine the two tests? I realize that a glm or other approach to analyze trends is the way to go, Any suggestions or insights would be appreciated.               ##### Attempt at using glm model ######         ltProb <- 0.4                       ## longterm average of nest survival probability         change <- c(.01,.03,.05)         yrs <- seq(1,20, by=1)                                                    ## Years of inquiry, probability of detecting a yearly change nest survival          samplesize <- 50                                               ## Reasonable sample size range         reps <- 1000                                                                       ## of simulations per          SurvProb <- ltProb                                                           ## initiating survival probablity for later use, nuisance         power <- matrix(nrow=length(change)*length(yrs)*length(samplesize), ncol=5) ## creating a matrix to hold data              scenario <- 0                                                                      ## initializing scenarios which will be a placeholder later on         for (a in 1:length(change)){                                                         ## loop through yearly pop change scenarios           for (c in 1:length(samplesize)){                                                                ## loop through sample size scenarios             for (b in 1:length(yrs)){                                                                              ## loop through years sceanarios               scenario=scenario+1               power[scenario,1] <- change[a]                                                                                ## filling matrix with pop change used               power[scenario,2]<-ltProb*((1-change[a])**yrs[b])               power[scenario,3] <- yrs[b]                                                                                ## filling matrix with number of years used               power[scenario,4] <- samplesize[c]                                                                                ## filling matrix with sample size used             }           }         }         colnames(power) <- c("PopChange", "Proportion2","yrs", "Sample_Size", "Power")         power=as.data.frame(power)         power$Sample_Size=as.numeric(power$Sample_Size)              scenario=levels(as.factor(power$PopChange)) ## To subset power             ps = levels(as.factor(power$Sample_Size))         results <- matrix(nrow=0, ncol=5) ## final matrix         results=as.data.frame(results)         reps=1000 ## set number of reps              for (k in 1:length(scenario)){           for (m in 1:length(ps)){             sub=power[(power$PopChange==scenario[k]) & power$Sample_Size==ps[m],] ## Subset by scenario and sample size             probs = matrix(nrow=1000,ncol=16) ## this is matrix for probabilities.             dat=matrix(nrow=0,ncol=3,NA)             dat=as.data.frame(dat)             for (l in 1:reps){ ##This should be for the replicates               dat=matrix(nrow=0,ncol=3,NA)               dat=as.data.frame(dat)               print(l)               for (j in 1:nrow(sub)){                 tmp=rbinom(n=sub$Sample_Size[j],size=1,prob=sub[j,2])                     tmp.dat=cbind(tmp,rep(sub$yrs[j],sub$Sample_Size[j]),rep(sub$Proportion2[j],sub$Sample_Size[j]))                     dat=rbind(dat,tmp.dat)                     if (sub$yrs[j]>4){                   x=(glm(dat[,1]~dat[,2],family=binomial))                   probs[l,j-4]=ifelse(summary(x)$coefficients[2,4]<.05,1,0)                     }                   }                 }                 sub$Power[5:20]=apply(probs,2,sum)/nrow(probs)             results=rbind(results,sub)           }         }              tmp=results[!is.na(results$Power),] ### remove NAs             tmp$SampleSize_Decline=paste(tmp$Sample_Size,tmp$PopChange,sep=":")         tmp$SampleSize_Decline=as.factor(tmp$SampleSize_Decline)              ### Now do the same using prop.test #####         tmp$PropPower = NA             for (i in 1:nrow(tmp)){               print(i)               x=0               for (j in 1:1000){                 a = rbinom(tmp$Sample_Size[i],1,prob = (tmp$Proportion2[i]))                 b=binom.test(sum(a),length(a),p=.4,conf.level=.95,alternative="less")                 if (b$p.value < .05) {x=x+1}               }               tmp$PropPower[i] = x/1000         }              ##### graph results         d=unique(tmp$PopChange)             for (i in 1:length(d)){               par(mfcol = c(1,2))               sub = tmp[tmp$PopChange == d[i],]               plot(sub$yrs,sub$Power,col="red", main = d[i], xlab = "Number of years",ylab="Power")               points(sub$yrs,sub$PropPower,col="black")               plot(sub$PropPower,sub$Power,xlab="Binomial Proportion Power", ylab="GLM Power",main=d[i])           abline(0,1)           readline("Press <return to continue")         }                   ##### End Code #####


title: How fitting a mixture distribution of noncentral Student t components to a one-dimensional sample in R?

text: I`d like to extract the parameters of a two-component mixture distribution of noncentral student t distributions which first has to be fitted to a one- dimensional sample. My question is closely related to this thread, but as pointed out I want to use Student t components for the mixture: Which R package to use to calculate component parameters for a mixture model There are many packages for R that are capable of handling mixture distributions in one way or another. Some in the context of a Bayesian framework requiring kernels. Some in a regression framework. Some in a nonparametric framework. ... In general the "mixdist"-package seems to come closest to my wish. This package fits parametric mixture distributions to a sample of data. Unfortunately it doesn`t support the student t distribution. I have also tried to manually set up a likelihood function as described here: http://stackoverflow.com/questions/6485597/r-how-to-fit-a-large-dataset-with- a-combination-of-distributions But my result is far from perfect. The "gamlss.mx"-package might be helping, but originally it seems to be set up for another context, i.e. regression. I tried to regress my data on a constant and then extract the parameters for the estimated mixture error distribution. Is this a valid approach? But with this approach the estimated parameters seem to be not directly accessable individually by some command (such as fit1$sigma). And more importantly there seem to be serious estimation problems even in pretty simple and nonambiguous cases. E.g. in example 2 (see syntax below) I simulated a mixture which looks like this: ![kernel density estimate of the mixture](http://i.stack.imgur.com/MG7AA.jpg) When trying to fit a two-component student t mixture to these data either I get this error message (the deeper meaning of which I don't understand): ![enter image description here](http://i.stack.imgur.com/UPvg4.jpg) or I get wrong results (convergenve is reached only after approximately two hours as can be seen from the output): ![enter image description here](http://i.stack.imgur.com/HjlfW.jpg) The means could be estimated well, but both the variance and the degrees of freedom are estimated badly. In the TF2 implementation of the student t, the sigma parameter denotes the standard deviation. Its estimate is NEGATIVE for the first component! And for the second component the degrees of freedom estimate is also NEGATIVE. Probably one should not use these results in practice :( By the way: Is there a way to restrict these degree-of-freedom coefficient estimates to be natural numbers? The following syntax is my gamlss.mx-setup so far:               library(gamlss.dist)     library(gamlss.mx)     library(MASS)          # example 1 (real data):     data(geyser)     plot(density(geyser$waiting) )     fit1 <- gamlssMX( waiting~1,data=geyser,family="TF2",K=2 )     fit1     # works fine          # example 2 (simulated data):     N <- 100000     components <- sample(1:2,prob=c(0.6,0.4),size=N,replace=TRUE)     mus <- c(3,-6)    # denotes the mean of component 1 and 2, respectively     sds <- c(1,9)     # ... the standard deviations     nus <- c(25,3)    # ... the degrees of freedom     mixsim <-data.frame(rTF2( N,mu=mus[components],sigma=sds[components],nu=nus[components] ))     colnames(mixsim) <- "MCsim"     plot(density(mixsim$MCsim) , xlim=c(-50,50))     fit2 <- gamlssMX(MCsim~1,data=mixsim,family="TF2",K=2)     fit2     # error message or strange results (this also happens when using a sample of S&P500 returns)      I would be very grateful for any advice! I've read through many related manuals and vignettes so far but I`m still lost. Thanks a lot in advance!! Jo


title: Is common factor analysis ever performed using the covariance matrix?

text: Common factor analysis entails the eigendecomposition and subsequent interpretation of $\mathbf{C}$ which is often described the correlation matrix with the diagonal elements replaced with approximate "communalities"; $\mathbf{C} = \left(\mathbf{R} - \text{diag}(\mathbf{R}^{+})^{+}\right)$, where $\mathbf{A}^{+}$ indicates the Moore-Penrose inverse of matrix $\mathbf{A}$, possibly with subsequent rotations. Is common factor analysis ever based on the eigendecomposition of a transformation of the covariance matrix, rather than the correlation matrix? (And I don't mean by first transforming the covariance matrix into the correlation matrix.)


title: How to write custom model in SAS Enterprise Miner using SAS Code Node?

text: I want to write a simple GAM model in SAS Enterpise Miner and then compare it with Neural Networks I made before using built-in methods. How can I do that? Of course I want to take advantage of the already partitioned set into learn, validate and test subsets. Let's assume I want to use the following code ( source: link ) :               PROC GAM data=kyphosis;         model kyphosis = spline(NumVert,df=3)                          spline(Age,df=3)                          spline(StartVert,df=3)                          /dist=logist;         output out=estimate p uclm lclm;     run;      I'd be grateful for an explaination how to fit that code into SAS Code node in SAS Enterprise Miner, so that I could access lift curve and compare this model's results with other ones.


title: Can anyone clarify the concept of a "sum of random variables"

text: In my probability class the terms "sums of random variables" is constantly used. However, I'm stuck on what exactly that means? Are we talking about the sum of a bunch of realizations from a random variable? If so, doesn't that add up to a single number? How does a sum of random variable realizations lead us to a distribution, or a cdf / pdf / function of any kind? And if it isn't random variable realizations, then what exactly is being added?


title: using t-test to compare multiple slopes?

text: I work with an escape response. Basically, I want to find differences between animals that should run and animals that stay still when a stimulus is applied to them. I've been trying several approaches related to velocity and normalization (due to non-normal data) but lately I've been thinking there could be an easier way to analyse the data. My crude data is position as a function of time. Here's a plot of 28 animals for a specific time lapse with linear regression for each case. ![enter image description here](http://i.stack.imgur.com/gqfEo.png) There's some animals that doesn't look nice adjusted to this lm but I think I can fix that by using total movement instead of position (if an animal is located at 3 and moves to 6 but then comes back I should have a total movement of 6 instead of zero). Since this is only one trial of stimulation and I have several ones to do I was wondering if it is correct to use the individual absolut slopes as the "only" value for each trial for each animal. If it is, then one could use a simple t-test find differences. Besides, I was trying to figure out a correct way to discard outliers (differentiate zero from noise and discard "big runners"). I think this approach might help but I don't know if it is correct. Here's the table with slopes, just in case. Maybe the linear model is not the best one and I should come back to nl-regressions and fighting with initial values.               variable slope value     UT1      -0.14     UT2      -25.59     UT3      8.36     UT4      62.73     UT5      -2.68     UT6       5.18     UT7      2.23     UT8     -18.41     UT9      10.91     UT10    -22.55     UT11     -5.73     UT12     -10.91     UT13       3.45     UT14      8.82     TR1       17.86     TR2      42.23     TR3      -2.00     TR4       6.55     TR5       1.36     TR6      4.82     TR7     -1.00     TR8      0.00     TR9      -58.00     TR10       5.77     TR11       0.00     TR12    -1.18     TR13    24.09     TR14    -0.55      Thank you very much


title: Generalized Linear Model instead of log transformation

text: I'm studying factors affecting the accuracy of taxpayers' perceptions of their own marginal income tax rates. The dependent variable (DV) is the misperception in percentage points. It ranges from 0 (i.e. accurate perception) to 35 (i.e. 35 percentage points mis-estimated). The frequency distribution is positively skewed because many more individuals make small estimation errors while fewer make larger errors. The independent variables (IV) are mainly categorical (e.g. whether the individual prepares his/her own tax return) except for one IV (taxable income) which is scale data. Initially, I had planned to use a log-linear model, with transformation of the DV by taking log (y+1) to deal with the zeros in the DV. However, I've read from many posts here and elsewhere that this is not ideal. Besides, the regression coefficients will not be amenable to an intuitive interpretation. Could anyone advise if this analysis could be carried out using a Generalized Linear Model instead or am I totally off-tangent here to consider that? Any pointers to set me in the right direction are very much welcomed!


title: Concerns related to neural network matlab toolbox

text: I have some concerns related to the use of `nntool` in MATLAB toolbox. Following links like this, I have found that `nntool` by default normalizes the inputs to the range `[-1 1]`. So I am a bit concerned, I created a neural network with `tansig` activation in the first layer and `logsig` activation in the output layer. I manually normalized the outputs to the range of `[0 1]` in the data and fed it to `nntool`. Now my question is does `nntool` further normalizes it to the range `[-1 1]`. If it does then it is not correct, the output of `logsig` cannot be in the range of `[-1 1]`. I am using the `newfit` function like this               load house_dataset;     net = newfit(houseInputs,houseTargets,20);      Any suggestions?


title: Ordering cross-sectional residuals for Durbin-Watson test

text: I'm working on a cross-sectional regression model. Is it correct to test the presence of specification error in the model using Durbin-Watson test (or the non-parametric Runs test) with residuals **ordered according fitted-values**?


title: How to handle differences of neocortical parcellations between groups

text: In my lab, we are currently looking at how two handedness groups (mixed handed and consistent handed) differ in terms of lGI (local Gyrification Index, a cortical folding measure) across 74 different brain regions. We are thinking of different methods of comparing the handedness groups. One idea that was proposed was to calculate mean lGI for each brain region (74 different means) for each group, and then compare across the groups. The proposed idea was to do a paired sample t-test between the handedness groups, the justification being that mean lGI for Region 1 for one hand group is paired to mean lGI for Region 1 for the other hand group. Two questions:   1. Is it justifiable to use a paired sample t-test in this context (and why or why not)?   2. Can you think of a way to approach this problem without  a. data reduction,   b. a multilevel model with 74 level 1 equations, or   c. 74 different regression equations? Note that we are also interested in seeing how sex and intracranial volume influence lGI.


title: Repeated Measures Factorial, I think - help?

text: I have complete data set that fits nicely into excel. One factor is rows (R). The other factor is columns (C). Another factor is related to the analysis of those (F3)... So my experimental design was to subset by row(R), column(C) and the other factor (F3 - has 2 levels). So all of my combinations of R,C and F3 are related to each other because they are subsamples from the same data. The subsample sizes of R and C get pretty close to their respective population sizes. So the higher levels of factor R and C are more related to each other than the lower levels. So the variance differs between levels based upon sample size and effects.   1. Can I analyze this with a repeated measures factorial?   2. How so? Can you point me in the right direction for some highly relevant reading? I haven't studied these types of designs before.


title: Does Newman's network modularity work for signed, weighted graphs?

text: The modularity of a graph is defined on its Wikipedia page. In a different post, somebody explained that modularity can easily be computed (and maximized) for weighted networks because the adjacency matrix $A_{ij}$ can as well contain valued ties. However, I would like to know whether this would also work with signed, valued edges, ranging, for instance, from -10 to +10. Can you provide an intuition, proof or reference on this issue? Thanks.


title: Test of significance of multiple proportions in two groups

text: I have a dataset of two genetically different cell types, "A" and "B", and each of them can have four different morphologies, let's say "spiky cell", "elongated cell", "round cell", "triangular cell". This yields a 2x4 contingency table:                             Spiky   Long    Round    Triangular     Type "A":       1      10      13       44         Type "B":      10       3       3       20         How do I test if the differences in morphologies are different? As in "Cell group B has a significantly higher proportion of spiky cells than cell group A". Normally, I would use a Fisher's exact test, and I have thought about just doing four Fisher's tests. However, since there are four different morphologies, I feel like I should be using a statistical model that models all four proportions. Like in the above example, it seems that there are more spiky "B" cells and more elongated "A" cells, and they are somehow not independent, since a cell can only be one of the four. Should I be using some kind of multilevel model? I have tried searching for answers in previous questions, but I didn't find a solution. I hope someone can suggest a good test or maybe direct me to a previous discussion. Thanks!


title: Interaction: ignore two-level if three-level is significant in every case?

text: I am using lmer in R to run LMM. My DV is continuous and my IVs are categorical. Many statistician said if the three-level parameter is significant, I cannot interpret the two-level parameter. Does it also include when the parameters are not contradict to each other too? For example, in the two-level, Group2:color2 is significant, but in three- level, Group2:color3:male is significant. Does it mean that I can interpret both levels as they are not contradict to each other, i.e. The first one is telling me that group2 and the other group is different in color2 than the baseline color, whereas the second one is telling me that the effect of group2 and color3 is higher in male compared to female. Or should in just ignore the two-level interaction regardless of whether than they contradict to each other or not?


title: stick breaking model of Dirichlet process

text: I have a question regarding sticking-breaking model of Dirichlet process, which is defined as follows: ![enter image description here](http://i.stack.imgur.com/8llR3.png) There are further statements that ![enter image description here](http://i.stack.imgur.com/9jKL5.png) I am not clear that how to derive equation 1 from that posterior distribution and why does the equation 1 is equal to equation 2. Here $\mathcal{B}$ is Beta distribution.


title: model compare by anova test

text: I did model test of in anova. once mod1 vs mod2 , other mod2 vs mod1. Does it matter? CAn I compare two models of different df and different parameters?               anova(reg5,reg4, test="Chisq")          Analysis of Variance Table     Model 1: Gas ~ CODload     Model 2: Gas ~ Flow * CODin       Res.Df   RSS Df Sum of Sq  Pr(>Chi)         1     65 27786                                2     63 12835  2     14952 < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          > anova(reg4,reg5, test="Chisq")          Analysis of Variance Table     Model 1: Gas ~ Flow * CODin     Model 2: Gas ~ CODload            Res.Df   RSS Df Sum of Sq  Pr(>Chi)         1     63 12835                                2     65 27786 -2    -14952 < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Comparison of Time Series Data

text: I'm new here, and fairly new to statistics, so I apologize if this is an easy question, but I'm stumped. I'm doing a research project simply for fun (yes, I'm weird like that). I was never able to find a professor at my graduate school interested in sponsoring the study so IRB refused to sign off on it so I was left to conduct it independently. I am somewhat stymied on expressing my finding mathematically and would greatly appreciate any help that you might be able to provide. My study analyzed 241,000 police records in two neighborhoods over a 17 year span. The crime data was compiled annually and plotted out. The first 12 years of the study one neighborhood had the variable that I was testing in effect, the last 5 years the variable was removed. This allowed me to have a baseline to compare the test neighborhood against. If the variable did not impact the crime rate, then I would expect to see that the two neighborhoods mirrored each other the entire time, consistently, which definitely happened in some cases. However, in other cases the result was a clear and drastic change between the control period and the test period, as pictured below. ![enter image description here](http://i.stack.imgur.com/zjEYt.png) I've tried doing trend line analysis, but I'm not sure how to express that difference mathematically. Another researcher suggested that I look at a t-score analysis, but the problem with that is that the t-scores only examine the mean values for comparison. In the case of this particular result pictured here, the actual mean value didn't change, but the year-to-year variance certainly did. Can anyone offer some advice on how I can be express this mathematically? Thanks in advance!


title: What is the sampling distribution of a sum of scores?

text: Can you help me understand, or refer me to a source about, what is the sampling distribution of the sum of scores? I know how to get the sampling distribution of the mean to construct a t-test, etc. More specifically how does the variability of the sampling distribution of sums of scores relates to the variability of the underlying variable in the population. e.g. People pass sand bags to each other. I measure the total time for a sample of n=10 people to pass the bag, and calculate an average. If I repeat this sampling I will have a standard deviation of the sample mean, and then I can calculate the variability in the population. Assuming a relatively normal process, even though I am measuring time, etc. The problem is if I measure the same sample of 10 people passing bags every time. The samples are related. Any advice?


title: Confusion related to histogram density estimation

text: I have some confusion related to how the density is estimated from the histogram. I have attached the screenshot of the paper as well. Any insights ![enter image description here](http://i.stack.imgur.com/AsOHR.png) I didn't get why you divide it into cubes and why is N=(1/h)^d and the formula of the density estimator how did it come?


title: Transforming a continuous covariate into a discrete one in logistic regression

text: I am largely inexperienced in the area of logistic regression. I was wondering if there were a good way to transform a continuous covariate in a logistic regression into a discrete one by subdividing the support of this continuous covariate in a smart way. Any reference is welcome if the question is too general.


title: Suggestion for statistical analysis. Is generalized estimating equation (GEE) a good option?

text: I have two groups a G1 and control. For both groups we measured their cell counts every 10 minutes: 0 to 650 min. We are trying to find:   1. Difference (in cell counts) between G1 and controls - significant   2. Effect of time (on the cell count) the same for both G1 and controls What statistic analysis do we need to use? One of my friends told me to use the generalized estimating equation (GEE), but I am not sure. What is GEE used for?


title: How to combine multiple similarity measures?

text: I have a hyperspectral image where the pixels are 21 channels. So each pixel $\in \mathbb{R}^{21}$. I want to perform clustering on the pixels with similarity defined by two different measures, one how close the pixels are, and the other how similar the pixel values are. Thus if $X_1$ and $X_2$ are the locations of pixels $p_1$ and $p_2$ I have: $$S_X = \|X_1-X_2\|^2_2$$ and $$S_p = \|p_1-p_2\|^2_2$$. I have seen these measures combined into a single measure like this: $$ S= e^{-\frac{S_p}{\sigma^2_p}} \times \,\, e^{-\frac{S_X}{\sigma^2_X}} $$ **My question:** Is there a right way and a wrong way to combine measures like this, or if it improves my clustering can I combine the measures in any way that suits me? My question is vaguely related to Combining multiple similarity measures.


title: How can I evaluate my model when the the testing data are too few and the generated results don't match testing data

text: I have a bioinformatic data set includes a very large negative examples (let's say 30000 examples) and just a few positive examples (let's say 150 positive examples). Since I need to feed an enough amount of positive examples (let's say 135 examples) to learn a model. I just have a few positive examples left for testing set. The problem is the positive predicted results don't match at all with the amount of positive examples (15 positive examples) in testing set. So are there some ways I can evaluate my learning method to know whether it's truly good or not? Thanks very much for your comments.


title: Precision - Recall: Graphical Representation

text: I'm a little bit confused with precision recall. I read some papers about recommender systems, where in one paper they have a graphical representation and in other papers they don't (they just have the values for precision/recall). What is the advantage of a graphical representation? Am I right, that for a graph you have to do something like a step function? e.g. you have a list of values `true,true,false,true,false` and you begin with the first one, calculate precision recall and then draw that point into the graph. For the next point, you take the first and the second value `true,true` into account, calculate precision recall, draw the point... and so on. Am I right that this is similar to drawing a ROC curve?


title: Rank items based on different parameters

text: Let's say I have **N** items to be ranked based on e.g. 5(it could be few or less) different performance measurements. And ranking must reflect combined effect of all perfornance parameters. Which stastical model should I use ? I know this is too little information. But, I don't know what else I should write. I can edit question as and when more clarfication asked by you guys. EDIT: Let me elobrate more and get into specifics. I want to rank my VoIP carriers in my application to optimize my call routing based on two performance parameters.               1) ACD(Avg Call Duration) = Total Duration / No of Successful calls     2) ASR(Avg Seziure Ratio) = No of Successful calls * 100 / Total Calls (Including failed calls)      I measure above parameters from time to time from accumulated past Call Detail Records. Usually the sample period is 6 hours. I may consider ranking individual dimensions and then averaging them. But,I don't want a particular dimension to dominate. That means, if a carrier provide good ACD but bad ASR - it's bad for me. Vice versa.


title: proposal distribution on a manifold

text: I try to sample a set of parameters $(x_i)$ from a posterior $p((x_i)|(o_j))$ using a metropolis-hastings algorithm. The fact is that a simple independent random walk on each of the $x_i$ as proposal leads to poor convergence results. To alleviate this, I would like to reparametrise the $(x_i)$ as: $(x^c_i),f$ with $\sum_i x^c_i=1$, in order to design a proposal of the form $$ p((x^c_i),f|(\tilde{x^c}_i),\tilde{f}) = p((x^c_i)|(\tilde{x^c}_i)) p(f|\tilde{f}) $$ for which I expect better convergence properties. My question concerns the way to design a $p((x^c_i)|(\tilde{x^c}_i))$ ensuring a tractable sampler; I did not manage to find any reference about that. Thanks. ps: I am not very confident with the title of the post


title: Why is a 'Markov Random Field' a field?

text: In particular, I'm not seeing how a commutative group is defined on a graph, nor addition and multiplication. So why is it a field?


title: How to check if my regression model is good

text: One way to find accuracy of the logistic regression model using 'glm' is to find AUC plot. How to check the same for regression model found with continuous response variable (family = 'gaussian')? What methods are used to check how well does my regression model fit the data?


title: Estimating age based on height

text: I wonder if it's possible to estimate a child's age given the child's measured height. I have found this height chart: http://resource.nhi.no/resource/4281-21-hoyde- gutter-5-19-who.pdf Is it possible to use the information in this chart to estimate what the probability for a child to be of a certain age when the height is known? I'm sorry, but my knowledge in statistics are very small.


title: How many attributes to select for classification

text: I am aware that this question is vague, but after reading multiple standard books of the community, I am still wondering about the following applied problem. Lets say I have a dichotomy problem and for training I have $M$ observations, each consisting of $k$ features (attributes). I would be interested in any literature which helps determining the number of features one should select for the classification process.


title: Dummycoding based on clustering from OM distances

text: I'm using TraMineR to determine a certain clustering based on Optimal Matching distances:               # Define a sequence object     sequences.seq <- seqdef(sequences, left="DEL", right="DEL", gaps="DEL", missing="")          # Compute OM costs     costs <- seqsubm(sequences.seq, method="TRATE")          # Create OM object     sequences.om <- seqdist(sequences.seq, method="OM", indel=1, sm=costs, with.missing=FALSE, norm="maxdist")          # Use OM object to create a dendrogram     clusterward <- agnes(sequences.om, diss = TRUE, method = "ward")     plot(clusterward, labels=colnames(sequences))      This gives me a plot of a dendrogram. Now I want to dummycode the cases in my dataset dependent on which cluster they fall in. How can I accomplish this?


title: Multinomial Conditional Distribution

text: I'm working on a problem and after an 1.5, I've not had much of a progress. I appreciate if you can guide me on how to solve it. Consider $$X = {X_{ij}: i = 1, \dots, I, j = 1, \dots J} \sim Multinomial (n,p)$$ Also consider $X_{.} = \{X_{i.} = \sum_j X_{ij}: i = 1, \dots, I\}$. How can I show that $X|X_{.}$ is product multinomial meaning that each of the row vectors are conditionally multinomial and are mutually independent. Thanks for your help.


title: SVM in R package e1071

text: I am trying to use SVM to make a prediction (True or False) on a dataset with many independent variables. I am wondering how I can identify the most useful variable in making the prediction. I understand that there is a function that allows me to calculate the "weights" of the variables in R package `e1071`.               w=t(model$coefs)%*% model$SV      But I am not sure how I can decide which are the important variables, based on these weights.


title: Need help about a longitudinal study with eye tracking data

text: I am currently conducting a longitudinal research on learner reading development in second language. Let me explain it briefly. I am using eye tracking procedures and trying to describe the development of reading skills of my participants. Saccade lenght (how fast they read), fixation durations (how long they look at the text) and regressions (how many times they return and read again) are my main variables which predicts their reading skill level. For example, less fixation and less regression means higher reading ability. So this is a study which will continue two semester (about 8 months) and in this period, I will conduct about 6 eye tracking sessions after classroom instruction. These 6 sessions include all the variables I mentioned above and which statistical tests should I use to describe the developments in their reading ability? Indeed this is not a pre test post test issue but rather a development in behavior. Can I apply hierarchical regression and how? What do you think? Thank you in advance.


title: How should I design an experiment with non-binary factors?

text: We are preparing an experiment that will run during the coming winter, where we assess certain energy performance indicators on a set of single-family houses in Europe. We are interested in finding out which factors have the greatest impact on these indicators. These factors will include the age of the building, its insulation level, its inertia, etc. I'm minded to treat these factors as, well, factors taking on a set of discrete values. For example, the age of the building will be coded as the decade in which it was built; the insulation level will be factored as poor, medium, good; etc. I know more or less how to design an (fractional) factorial design when the factors can take on 2 values. But how should I design an experiment where, for the sake of the example, one factor can have 5 levels, another 3, and another 3?


title: Beta-t-EGARCH model, what is the unconditional variance

text: I am trying to implement Beta-t-EGARCH model in R, however, I fail to understand what is the unconditional variance of this model. I know for simple GARCH model, unconditional variance is               ARCH0/ (1-ARCH1 - Beta1),      but what I have read for Beta-t-EGARCH model the unconditional variance is equal to Omega / (1-phi1). So in my case it would be:               Coefficients:                       omega         phi1      kappa1   kappastar        df              Estimate:  -2.9495492 0.9990466008 0.040487373 0.007719683 8.5488304          -2.84/(1-0.99904) = uncond. variance      and as you can see it does not make sense. If anyone knows how to calculate unconditional variance, please help. Here is a link to a pdf with the model and R code


title: Solving linear regression with weights and constraints

text: I would like to solve a linear regression (in R) with weights $w$ and a constraint. In other words, I would like to find $x$ that minimizes the sum of squares $$\sum_i w_i(b_i-Ax_i)^2$$ On top of that I have an external vector $d$, which I would like to use in a constraint, such that $d \cdot x \le 5$. Is this something that would be possible to do in R with `solve.QP` or perhaps some other R script? **Edit** : I am adding a bounty for a solution that doesn't require any other custom software except the cran packages. While rstan works perfectly unfortunately I am unable to install it on my production servers due to old versions of some libraries.


title: Randomized Response With Continuous Outcome

text: Can someone point me to an example of how to do randomized response with a _continuous_ outcome $Y \in [ 0,\infty ]$, rather than a binary one? I am interested in the mean and variance of this outcome or the entire distribution if possible. The true distribution is expected to have a mass at zero and a long right tail, similar to expenditures on hospitalization in the US. Poole (1974)* developed an approach where the interviewee multiplies his true response $Y$ by a random number $X$ and only tells the interviewer the product $Z=X \cdot Y$. Unfortunately, I am finding the notation and the example very hard to follow. * * * *The reference in the link above is gated, but the full citation is:  Poole, W. Kenneth (1974) "Estimation of the Distribution Function of a Continuous Type Random Variable through Randomized Response," _Journal of the American Statistical Association_ , Vol. 69, Iss. 348, pp. 1002-1005.


title: Regularization in multinomial manifold (k-simplex) space

text: I have a multi-class classification problem. In order to reduce the feature space, I have mapped my weights ($w$) to a multinomial manifold (or k-simplex) and compute the new weights ($v$) as follows: $$ v_i = \frac{w_i}{\sum\limits_{i=1}^n w_i} $$ Clearly, in the new manifold space, my weights lie on a hyperplane. My question is how to do regularization in this new space? What would be the equivalent of $l_1$, $l_2$ regularizers here? Any suggestions would be helpful to me.


title: Understanding the Chi-squared test and the Chi-squared distribution

text: I am trying to understand the logic behind chi-squared test. The Chi-squared test is $\chi ^2 = \sum \frac{(obs-exp)^2}{exp}$. $\chi ^2$ is then compared to a Chi-squared distribution to find out a p.value in order to reject or not the null hypothesis. $H_0$: the observations come from the distribution we used to created our expected values. For example, we could test if the probability of obtaining `head` is given by $p$ as we expect. So we flip 100 times and find $n_H$ `Heads` and $1-n_H$ `tails`. We want to compare our finding to what is expected ($100 \cdot p$). We could as well use a binomial distribution but it is not the point of the question… The question is: **Can you please explain why, under the null hypothesis, $\sum \frac{(obs- exp)^2}{exp}$ follows a chi-squared distribution?** All I know about the Chi-squared distribution is that the chi-squared distribution of degree $k$ is the sum of $k$ squared standard normal distribution.


title: Understanding Conjugate families

text: I'm trying to understand conjugate families of distributions, and am hoping someone could give me a hint or direct to a good paper (or book) covering basics. Or perhaps, simple example would help?   Let's assume Gaussian model with both mean and variance random.   Uniform sample of the universe, say, is {2, 3, 4, 4, 5, 6, 8}, and sample of some subset is {4, 5, 5, 7}. (or any other values).   **1)** How hyperparameters should be estimated? I know that if variance is known, mean of the subset could be found as $\mu' = \frac{\tau_0 \mu_0 + \tau \sum_{i=1}^n x_i}{\tau_0 + n \tau}$. **2)** But what about estimating mean in our case?   **3)** And what about variance of the subset? Thanks in advance!


title: How to compare time series with cyclical data, and describe any changes or trends

text: I have a bunch of time series where the data has a natural (known) cycle, for example daily or annual (or both). Here is an example (this is 6 years worth of temperature data sampled hourly): ![enter image description here](http://i.stack.imgur.com/MzxpY.png) Source data: http://pastebin.com/GESpGqVu How can this data be compared? How can we check (for example) whether the first three years are significantly different from the last three years? That a particular month is significantly different from that month in all other years? That a particular day is unusual? (compared to both the same day in other years and adjacent days in the same year) I would like a test for change vs no change, something that works as the cyclical data equivalent to a test for stationarity (unit root test, Dickey–Fuller, KPSS, etc). I would also like some kind of descriptive statistics that works in this case. It is good to know that there is some change, but what is it? :) I am mainly interested in detecting trends (eg: the summer-winter swing is increasing over time) and outliers (eg: a particular month is unusually cold) while taking all the available information about cyclical variation into account. Some of the answers to other questions about time series comparison suggest doing ARMA/ARIMA etc on the two series (or parts of a series) that one is interested in comparing, and then an F test to see if the model parameters change. This may address the first part of my question (have the series changed), but it does not address the second part (how to explain in human- understandable terms how they have changed). I am also interested in any other (not autoregressive model comparison) standard or particularly good approaches to this. See also (related questions, no full answer):   * Comparison of Time Series Data   * How to compare two time series?   * How to statistically compare two time series?


title: GLM with multiple categorical variables in R : how to interpret the result?

text: I have a binomial variable that I regress against different categorical variables which I have contrasted to build a reference of an individual Female, Married, aged 35-45, High education :               Call:     glm(formula = wpti ~ gender.f + is.married.f + age.f + edu.f,      family = binomial(link = "logit"))          Coefficients:                        Estimate Std. Error z value Pr(>|z|)         (Intercept)        -1.303107   0.046522 -28.010  < 2e-16 ***     gender.fMale       -0.537958   0.036833 -14.605  < 2e-16 ***     is.married.fSingle  0.012196   0.040584   0.301 0.763792         age.f<25           -0.298081   0.078040  -3.820 0.000134 ***     age.f25-35         -0.121670   0.051283  -2.373 0.017667 *       age.f45-55         -0.006033   0.049078  -0.123 0.902162         age.f>55            0.239855   0.058727   4.084 4.42e-05 ***     edu.fLow           -0.340115   0.049879  -6.819 9.18e-12 ***     edu.fMedium        -0.298925   0.041729  -7.163 7.86e-13 ***      I was expecting to see all the coefficients (and not only the one for gender) to change if I contrast the gender factor with Female as reference, but they all remain the same, as if being a man or woman has no effect.               Call:     glm(formula = wpti ~ gender.f + is.married.f + age.f + edu.f,      family = binomial(link = "logit"))           Coefficients:                     Estimate Std. Error z value Pr(>|z|)         (Intercept)        -1.841064   0.048257 -38.151  < 2e-16 ***     gender.fFemale      0.537958   0.036833  14.605  < 2e-16 ***     is.married.fSingle  0.012196   0.040584   0.301 0.763792         age.f<25           -0.298081   0.078040  -3.820 0.000134 ***     age.f25-35         -0.121670   0.051283  -2.373 0.017667 *       age.f45-55         -0.006033   0.049078  -0.123 0.902162         age.f>55            0.239855   0.058727   4.084 4.42e-05 ***     edu.fLow           -0.340115   0.049879  -6.819 9.18e-12 ***     edu.fMedium        -0.298925   0.041729  -7.163 7.86e-13 ***      However, if I remove the gender in the regression, the coefficients for the other variables are different :               Call:     glm(formula = wpti ~ is.married.f + age.f + edu.f, family = binomial(link = "logit"))          Coefficients:                         Estimate Std. Error z value Pr(>|z|)         (Intercept)        -1.562816   0.043369 -36.035  < 2e-16 ***     is.married.fSingle  0.071780   0.040378   1.778 0.075455 .       age.f<25           -0.330773   0.077982  -4.242 2.22e-05 ***     age.f25-35         -0.134088   0.051091  -2.625 0.008677 **      age.f45-55          0.004803   0.048860   0.098 0.921694         age.f>55            0.224112   0.058419   3.836 0.000125 ***     edu.fLow           -0.411232   0.049424  -8.321  < 2e-16 ***     edu.fMedium        -0.332353   0.041488  -8.011 1.14e-15 ***      So I am confused : I was under the impression that the coefficients could show the difference of likeliness of having a result positive or negative if only ONE variable would vary. Now I dont know what to think. Can somebody explain how to interpret the results, and why it is different if I remove one categorical variable (eg. gender here) althought it seems it has no effect on the others ?


title: Estimation in Naive Bayes

text: I have a very silly question. In Multinomial Naive Bayes Classifier, which parameter estimation do we use, is it Maximum Likelihood or Maximum A Posteriori? If any one of the esteemed members may kindly help me out.


title: When can a pilot sample be included in the main study?

text: For example, when estimating the proportion of database records that agree with paper documentation. Auditing each record is tedious, so it is tempting to include the pilot sample in the main study (once the required main sample size is estimated). Provided the sampling method is the same, does this introduce bias?


title: Frequent Occurrence of Rotated Graph of Logistic Function

text: I have found that in much of the data that I am looking through, after sorting from largest to smallest, there is a pattern similar to a rotated logistic function. That is, it declines steeply along the vertical axis at first, evens out, and then continues to decline along the vertical axis at the same rate as it declined to start with. Does this tell me anything about the nature of the data underlying the graph outside the obvious? I feel like I might be missing something from this. If this is a ridiculous thought just let me know and I can remove the question.


title: how does Cox proportional hazards model deal with time-dependent variables

text: When considering time dependent data in survival analysis, you have multiple start-stop times for an individual subject with measurements for the covariates. If each season has a different size (for example: repr=90 days,post_r=5,winter=23), the probability of an individual dying in repr it's largest. How does the Cox model deal with different sizes of time intervals? I'm using `coxph()` in R. Here's an example:               subject | start   | stop   | event  | season    |       --------+---------+--------+--------+-----------+       1       | 1       | 90     | 0      | repr      |       1       | 90      | 95     | 0      | post-r    |       1       | 95      | 118    | 1      | winter    |       2       | 1       | 23     | 0      | winter    |       2       | 23      | 113    | 0      | repr      |       2       | 113     | 118    | 0      | post-r    |       2       | 118     | 141    | 1      | winter    |


title: Goodness of fit for generalized linear nodel

text: How to measure the "goodness of fit" in generalized linear models for repeated measures. To make clear: I'm searching for something which is based on the (deviance) residuals, something like an adj. R^2 for linear models. Something that could make me state about the efficacy of the model? Any thought would be appreciated!


title: Do reliability and construct validity have same technical(statistical) assumptions?

text: It appears that reliability and construct validity have a similar backing in terms of statistics. Is it true ?


title: What is the correct definition for completness in d-separation in directed graphical models?

text: I was reading Koller's book of probabilistic graphical models and in section 3.3.2 she discusses what properties should hold for d-separation as a method for determining independence. She tries to motivate that hopefully, separation should be sound and complete. _Sound_ in the sense that, if we have that two nodes $X$ and $Y$ are d-separated given some set of nodes in the graph $\vec{Z}$, then we are guaranteed that they are in fact conditionally independent given $\vec{Z}$. _Complete_ in the sense that hopefully, d-separation detects all possible independencies. She tries to make those terms rigorous by providing some explanations and definitions and I wish to ask about the definitions she provides for _completeness_. First define that $I(P)$ is the set of all independencies that hold in distribution P and similarly, $I(G)$ are all independencies that hold in graph G. Also, let the notation $d-sep_G(X;Y \perp \vec{Z})$ be defined as; X is d-separated from Y in graph $G$ given (the set of nodes) $\vec{Z}$. In that section (3.3.2) she defines the following notion: > A distribution P is **faithful** to G if, whenever $(X \perp Y | > \vec{Z}^{\,}) \in I(P)$, then $d-sep_G(X;Y \perp \vec{Z})$. In other word, > any independence in P is reflected in the d-separation properties of the > graph. The following is the important part. In trying to define for completeness, she suggests the following notion of completeness: > 1) For any distribution P that factorizes over $G$, we have that P is > **faithful** ; if X and Y are not d-separated given $\vec{Z}$ in G, then X > and Y are dependent in all distribution P that factorize over G. However, it ends up that this is actually not a good definition for completeness (for some reason which I don't completely understand yet but is not the main issue in this question) and decides to define a **weaker** definition. The definition is instead: > 2) If $(X \perp Y \mid \vec{Z})$ in all distributions P that factorize over > $G$, then $d-sep_G(X;Y \mid \vec{Z})$. And the contrapositive: if $X$ and > $Y$ are not d-separated given $\vec{Z}$ in G, then X and Y are dependent in > some distribution P that factorizes over G So this is the main issue with my question. Even though she provides two supposedly different definitions, I fail to appreciate what is the true **difference** between the two. It seems to me that they both say the exact same thing (I know, they don't because she says that the second is "weaker"), but I wanted to understand the difference between the two definition and why the second one is weaker.


title: Measure score change over time while accounting for baseline differences

text: I'd like to test for and estimate group differences in NIHSS (National Institute of Health Stroke Scale) change between hospital discharge and three months after hospital discharge. Because the score is bounded between 0 and 42, the signed difference between the two measurements seems to me inappropriate, as people with a higher score at the first measurement cannot possibly improve as much as people with a worse score at that time, and vice versa. How can I account for this? I considered a logistic regression model with the change in scores as dependent variable (as the score can be interpreted as a proportion, ie. x/42), but it I don't think I can control in this model for the baseline score, as it is already contained in the change in scores. It is also of importance that the NIHSS is somewhere between interval and ordinal, being a compound of several different subscores. It is not linear. I have found some more or less related threads here (such as Should the difference between control and treatment be modelled explicitly or implicitly? or Is it valid to include a baseline measure as control variable when testing the effect of an independent variable on change scores?), and am aware of the 1990 paper by Paul Allison ( http://www.statisticalhorizons.com/wp- content/uploads/Allison.SM90.pdf ). Edit: What I have looked into


title: Testing significance of individual coefficients in multi-linear regression

text: Traditionally, the individual coefficient significance test can be solved by applying t-test. However, in my case, as follows:   1. the predictor values have been standardized   2. the standardized coefficient have been known   3. we would not know the values of predictors, thus we can not use the t-test I argue that whether i can directly use the standardized coefficient to derive which one is significant, such as using confidence interval. Thanks!


title: What is a reasonable independence test for a time series?

text: Suppose that I have a time series, $x_{1:n}$, and I would like to test whether they are all independent samples from the uniform distribution on $[0,1]$. This is a somewhat basic question, but how do I correctly test for that? I am using the Kolmogorov-Smirnov test to check that they are uniform. One thing I can do is calculate the autocorrelation function for $x_{1:n}-\frac12$, which should be small. It should lie within approximately $\pm 2/\sqrt{n}$ and I can calculate the Ljung-Box and Box-Pierce statistics. My only concern with them is that if I repeatedly calculate them for just known-random data, e.g., `rand(1000,1)`, the distribution of p-values for both of them is not quite uniform. Which, if I understand correctly, means they are not very good tests. What would be a good test for independence?


title: Sampling distribution of average of some covariance matrices

text: I have $K$ datasets, each with $N$ variables and $M$ samples (they are in fact EEG time series, but I discard time and treat them as $K$ iid multivariate samples) and assume they are coming from the same multivariate normal distribution. I am interested in estimating the covariance matrix. Now it can be done in two ways:   * Concatenating the datasets together, and calculating the covariance matrix. Its sampling distribution would be Wishart, given the assumptions of multinormality of samples.   * Calculating the covariances separately for each dataset and averaging those matrices (with arithmetic mean) to form one total covariance. The first method is straightforward and have well established properties, but the second is in most cases much more feasible in my environment. From properties of variance of Wishart distribution $\Big($S.E. of element $C_{i\,j}$ of covariance matrix $C$ equals $\sqrt{\frac{C_{i\,j}^2 + C_{i\,i} C_{j\,j}}{M-1}}$ $\Big)$ and CLT (Central Limit Theorem) I can see, that both expected value and standard error of estimate of covariance matrix should agree for both methods. But yet, (obviously) the methods don't generate numerically the same covariances.   1. Is it really true, that both ways of estimating covariance matrix have the same standard error?   2. Does Wishart distribution behavior can be approximated by normal distribution when sample size parameter goes to infinity (just like we do for chi squared distribution)? If so, what are the conditions to have reasonably good approximation?   3. I bet, that if one can approximate matrix elements of Wishart distribution with normal distribution, than the validity of the second method depends on validity of this approximation. But please, can someone correct me, if I am wrong? I need those answers to justify (if it is justifiable at all :-) ) interchangeable use of both estimators in the article about performance of some joint diagonalization of $C$ algorithms.


title: Comparison between MAD and SD

text: I am reading Huber's Robust Statistics (2nd). On page 2 and 3 he gave an example. The basic facts are summarized here. Let $(X_n)$ be a sequence of random variables and define two measures of spread as follows.   1. Mean Absolute Deviation: $d_n := \frac{1}{n}\sum|x_i-\bar x|$.   2. Standard Deviation: $s_n := \sqrt{\frac{1}{n}\sum (x_i-\bar x)^2}$. Then he mentioned that Fisher claimed that for identically distributed normal observations $s_n$ is about 12% more efficient than $d_n$. In addition, $s_n$ converges to $\sigma$ while $d_n$ converges to $\sigma\sqrt{2/\pi}\doteq 0.8\sigma$. I have several questions about these statements.   1. How to prove that $s_n$ is **12% more efficient** , please? As least where to find the proof, please?   2. How to prove that $d_n$ converges to **$\sigma\sqrt{2/\pi}\doteq 0.8\sigma$** , please? Again at least where to find the proof, please?   3. I did some simulation to test all the above statements. Here are the codes and outcome. > >     n <- 10000 # number of samples >     x <- array(list(), n) >   >     set.seed(2014) >   >     for(i in 1:n){ >       x[[i]] <- rnorm(10000) # the 10000 here is the size of each sample >     } >   >     dn <- rep(0, n) # mad >     sn <- rep(0, n) # sd >   >     for(i in 1:1000){ >       dn[i] <- mean(abs(x[[i]]-mean(x[[i]]))) # mad >       sn[i] <- sqrt(var(x[[i]])*999/1000) # sd >     } >   >     mean(dn) # **0.07979068 check out** >     mean(sn) # **0.09995901 check out** >   >     var(dn)/var(sn) # **0.6371817** >   As the above simulation shows, the 12% efficiency of $s_n$ does not check out. Why is this the case, please? Did I make errors in my simulation, please? Thank you!


title: Why did Thomas Bayes find Bayes' theorem so challenging?

text: This is more of a history of science question, but I hope it's on-topic here. I've read that Thomas Bayes only managed to discover Bayes' theorem for the special case of a uniform prior, and even then he struggled with it, apparently. Considering how trivial the general Bayes' theorem is in modern treatment, why did it present a challenge for Bayes and other mathematicians at the time? For comparison, Isaac Newton's Philosophiæ Naturalis Principia Mathematica was published 36 years before Bayes' main work, which was published posthumously.


title: Stores and ordered rank of products

text: I have data about popularity of some products in some stores. For each of the stores I have only a ordered rank of best-selling products in them. I need help how to find information about which products are popular only in certain stores versus products that are popular in all stores.


title: Simple question on weak law of large numbers

text: ![enter image description here](http://i.stack.imgur.com/1amt2.png) What I have so far is this: I need to find the E[.] and V[.] and show that the expected value equals to $\sigma^2$ and variance tends to 0 as T -> $\infty$. I already found the Expected value to be equal to $\sigma^2$. How do i show the variance tends to 0? Is my approach right?


title: What is a calibrated statistical test

text: I am a little bit confused about the concept of "well calibrated **tests** " -- what aspects are we referring to from **a test** if it is well calibrated? I tried to search some papers on such discussions or textbooks, but the results seem to be very subject-specific. Thank you for your clarifications.


title: Descriptive and Inferential vs Parametric and Non-Parametric Statistics

text: Are the concepts of descriptive vs inferential statistics and parametric vs non-parametric statistics orthogonal? As in, can we have a descriptive parametric statistic or a descriptive non-parametric statistic? if so what are some examples of? I know you can have an inferential parametric or non- parametric statistic. My original thought was that parametric vs non-parametric statistics falls in a subcategory of inferential statistics, and that inferential statistics makes use of descriptive statistics. but some sources I read seem to say that you can have descriptive non-parametric stats. Anyone can shed some light on this?


title: Modeling a list with a tunable degree of disorder/shuffling

text: Imagine we have a list of ordered numbers $L = (1, 2,\dots, N)$. I want to add an arbitrary amount of "disorder" to that list. For instance:   * Adding a little bit of disorder would permute a few neighboring numbers   * Adding more disorder would permute several numbers, many of them not neighbors   * Adding maximal disorder would permute all the numbers Ideally I would like to use an algorithm that has a sound statistical/information theoretic foundation. I am not sure about how to search for this. I have tried Googling things like:   * Tunable degree of shuffling   * Generating random lists of numbers with a given entropy / Kullback-Leibler distance   * Generating a list with a given level of presortedness Note that the solution is not just shuffling a given percentage of the numbers. When disorder is small, the shuffling should occur among numbers that are close. Any ideas? Many thanks!


title: What method to use to analyse 3 groups of independent variables?

text: I’m trying to find a better way to analyse my data. Unfortunately my limited background in math does not allow me to freely make a confident choice. I’ve already done an ANOVA, the method mostly used in our lab. As it often happens with statistics in non-statistical fields, it was chosen mainly based on the criterion, “Everybody uses it,” rather than, “It is the best way to do it.” This approach is leaving me in a state of an itch. I understand, ANOVA indeed could be the best choice in my analysis, but again, I am not capable enough to make comparisons between different methods. I would appreciate the opinions and advice of a statistic-minded community. * * * Short overview of my project. I’ll try to make it as less specific as possible. Let’s imagine a brick wall. It has standard bricks and mortar in it. The wall has a measurable quality, let’s say a ratio between how long it needs to be vibrated until it collapses and the amplitude of the vibration. Call it G. G depends on both wall components—bricks and mortar. Mortar has two components on its own, like cement and sand. What I am trying to find is what relationship my special shaped bricks may play on overall G in different conditions. * * * # Independent variables. I divided my independent variables into 3 groups—A, B & C. **Group A** Changes in the mortar composition and has following arrangements: A0—nothing changed A1—reduced cement from the mortar A2—increased cement in the mortar A3—cement replaced to the new type of the cement known to affect G of the wall A4—changed shape of the sand particles **Group B** Changes in the shape of the bricks B0—normal wall with all normal bricks B1—added bricks with Shape 1 into wall B2—added bricks with Shape 2 into wall B3—added bricks with Shape 1 and 2 together **Group C** Layer of paint on the wall C0—no paint C1—thick layer C2—thin layer * * * For every combination of A, B and C I have 20 to 50 readouts of G. What would the best method to employ be to say, e.g., "Combination of two shapes of the bricks together is the main factor in keeping G at a certain level while mortar has changed and paint is applied?" I was thinking to make a multiple regression analysis, as I also have quantitative data of A0, A1 and A2, but not A3 and A4. Making multiple ANOVA comparisons in this case feels the same as making multiple $t$-tests. I would appreciate any thoughts. * * * Edit Here is how my creative graphing looks like.![Effect B](http://i.stack.imgur.com/2yyPU.jpg) Orange belt in the middle is a “zone of insignificance” from my ANOVA (compare mean of each column with mean of control column (A0B0C0). What I wanted to say is something along this line: “The effect of C2 on the wall is attenuated while B3 is present despite disturbance of the mortar with A1-A4"


title: Why would all the tests for normality reject the null hypothesis?

text: The Kolgomorov-Smirnov test, Shapiro test, etc.... all reject the hypothesis that a distribution is normal. Yet when I plot the normal quantiles and and histogram, the data is clearly normal. Maybe because the power of the tests are high? The sample size is around 650. So shouldn't at least one of these tests fail to reject the null hypothesis? Results:                          Kolmogorov-Smirnov    D          0.05031          Pr > D       <0.010                Cramer-von Mises      W-Sq       0.30003          Pr > W-Sq    <0.005                Anderson-Darling      A-Sq       1.66965          Pr > A-Sq    <0.005                Chi-Square            Chi-Sq  3250.43596     18   Pr > Chi-Sq  <0.001


title: Date/Time data on the y-axis

text: I have some data values of type date/time (the last date that a resource was accessed) and I wish to chart this data on the y-axis against the different categories of resource on the x-axis. What would be a sensible type of chart for displaying this sort of data ? For example is a histogram / bar chart satisfactory or is there a different graph type more applicable ? EDIT: Are bar charts for only displaying quantities - can dates be considered quantities ?


title: Standard reference for classical mathematical statistics?

text: Can anyone recommend some books that are considered to be standard references for classical (frequentist) statistics? IE, fairly comprehensive, and also, been around for a while so that typos and mistakes in formulas had a chance to be checked and corrected


title: How should one approch Project Euler problem 213

text: I would like to solve Project Euler 213 but don't know where to start because I'm a layperson in the field of Statistics, notice that an accurate answer is required so the Monte Carlo method won't work. Could you recommend some statistics topics for me to read on? Please do not post the solution here.


title: What is the expected gain from taking a winner's bet?

text: Suppose a gambler makes a series of independent bets. In each bet, there's a 50% chance that he loses $\$100$. If he wins, he has an option to take $\$100$ or run a "winner's bet", instead of taking that reward. In the winner's bet, there is a 60% chance to win 50 dollars, a 30% chance to win 150 dollars, and a 10% to win 200 dollars. What should he do in order to maximize the expected reward? Note: This is not homework. I've been out of school for a decade. :)


title: Difference between ~lp() or simply ~ in R's locfit

text: (Please note the cross-post at http://stackoverflow.com/questions/7626347/difference-between-lp-or-simply-in- rs-locfit) I am not sure I see the difference between different examples for local logistic regression in the documentation of the gold standard locfit package for R: http://cran.r-project.org/web/packages/locfit/locfit.pdf I get starkingly different results with               fit2<-scb(closed_rule ~ lp(bl),deg=1,xlim=c(0,1),ev=lfgrid(100), family='binomial',alpha=cbind(0,0.3),kern="parm")      from               fit2<-scb(closed_rule ~ bl,deg=1,xlim=c(0,1),ev=lfgrid(100), family='binomial',alpha=cbind(0,0.3),kern="parm")      . What is the nature of the difference? Maybe that can help me phrase which I wanted. I had in mind an index linear in bl within a logistic link function predicting the probability of closed_rule. The documentation of lp says that it fits a local polynomial -- which is great, but I thought that would happen even if I leave it out. And in any case, the documentation has examples for "local logistic regression" either way...


title: Help interpret cv.glmnet behavior for very large set of potential predictors

text: I've created this toy example to demonstrate something that is occurring with my real data. I do not understand how to interpret the apparent failure of cv.glmnet to find a solution when it's explicitly "doped" into the data. Let's create a matrix of random predictors: xtrain = replicate(100,rnorm(10)) Now let's create a response variable that is explicitly a function of some of the x: ytrain = 1 + 2*x[,1] - 3*x[,2] Run cv.glmnet: cvfit=cv.glmnet(x,y) If we query the cvfit object, we'll find that lambda.min is not equal to lambda.1se, as expected, nor is lambda.min the largest magnitude lambda value. If we make predictions on the original training data ypred=predict(cv.glmnet,xnew=xtrain,s="lambda.min") we'll get quite good results in comparing ypred to ytrain. HOWEVER, if we do the exact same thing with a much larger set of potential predictors. xtrain=replicate(100000,rnorm(10)) ytrain=1 + 2*x[,1] - 3*x[,3] cvfit=cv.glmnet(xtrain,ytrain) we end up with a cvfit object where lambda.min = lambda.1se and this lambda value is the largest magnitude lambda value. The number of non-zero coefficients at that lambda is zero. And, of course, predictions at this lambda are nothing but the mean of ytrain. I could use some help understanding this limitation in applying the cv.glmnet to such a problem and if there is a work-around of some sort (e.g. taking random subsets of the xtrain matrix).


title: Can mean values be seen as a form of PCA with a LOT of assumptions?

text: This question just occurred to me out of the blue. PCA is a way to reduce dimensions. Another way that is often (perhaps too often) used is to take mean values of two or more variables. This is done a great deal to create "scales" or "ratings" that get published in the popular press about various political issues. Dimensions are naively "averaged" and this cumulative "average" is portrayed as a magical one-number-sums-all figure. Is this, nevertheless, a very crude relative of PCA? If so, what assumptions does it have? I can think of the following: Sphericity is perfect or near-perfect. All variables are perfectly orthogonal from all other variables. Variables are all perfectly scaled against each other. What else would there be?


title: Trying a multivariate analyses on time series (with R)

text: I got measures of one variable (that behaves as a time series) for different conditions (some quantitatives, but mostly are qualitatives). For example, this is a "fake" representative plot of this variable for 2 differents conditions : http://www.casimages.com/img.php?i=140415020040937841.png How can I statistically point out difference for my measured variable using Correspondence analysis methods ? Thaks for help


title: How to compute the confidence interval of the ratio of two normal means

text: I want to derive the limits for the $100(1-\alpha)\%$ confidence interval for the ratio of two means.   Suppose, $X_1 \sim N(\theta_1, \sigma^2)$ and $X_2 \sim N(\theta_2, \sigma^2)$ being independent, the mean ratio $\Gamma = \theta_1/\theta_2$. I tried to solve: $$\text{Pr}(-z(\alpha/2)) \leq X_1 - \Gamma X_2 / \sigma \sqrt {1 + \gamma^2} \leq z(\alpha/2)) = 1 - \alpha$$ but that equation couldn't be solved for many cases (no roots). Am I doing something wrong? Is there a better approach? Thanks


title: How to test homogeneity of variance if the data may not be normally distributed?

text: I need a test to check for homogeneity of variances. I am not sure if a series is normally distributed.   * If I use a test of homogeneity of variance that works for NON-normally distributed data, does the series also pass (with "pass" I mean that the variances ARE homogeneous) if the distribution is normal ?   * Can we say that a test of "non-normal" is more tolerant than a "normal" test? (obviously I'm always referring to the distribution of the data) I know there are tests to check if the distribution is normal or not, but I would like to take away this doubt.


title: Free econometrics textbooks

text: Are there any free econometrics textbooks available online for individual use? I'm aware of Bruce Hansen's first year PhD Econometrics textbook, but I'd be interested to hear if there are any other such resources. Note: This question is similar to an existing CV question, but I'm specifically asking for references that would be useful to econometricians and students of economics.


title: How to interpret GARCH parameters?

text: I use a standard GARCH model: \begin{align} r_t=\sigma_t\epsilon_t\\ \sigma^2_t=\gamma_0 + \gamma_1 a_{t-1}^2 + \delta_1 \sigma^2_{t-1} \end{align} and $a_t=\sigma_t\,\epsilon_t$ I have different estimates of the coefficients and I need to interpret them. Therefore I am wondering about a nice interpretation, so what does $\gamma_0$,$\gamma_1$ and $\delta_1$ represent? I see that $\gamma_0$ is something like a constant part. So it represents kind of an "ambient volatility". The $\gamma_1$ represents the adjustment to past shocks. Also, the $\delta_1$ is not very intuitively for me: It represents the adjustment to pas volatility. But I would like to have a better and more comprehensive interpretation of these parameters. So can anyone give me a good explanation of what those parameters represent and how a change in the parameters could be explained (so what does it mean if e.g. the $\gamma_1$ increases?). Also, I looked it up in several books (e.g. in Tsay), but I could not find good information, so any literature recommendation about the interpretation of these parameters would be appreciated. Edit: I would be also interested in how to interpret the persistence. So what is exactly persistence? In some books I read, that the persistence of a GARCH(1,1) is $\gamma_1+\delta_1$, but e.g. in the book by Carol Alexander on page 283 he talks about only the $\beta$ parameter (my $\delta_1$) being the persistence parameter. So is there a difference between persistence in volatility ($\sigma_t$) and persistence in shocks ($a_t$)? ![vo](http://i.stack.imgur.com/2t9hd.png)


title: Comparing densities of two partitions of a population

text: I have a population with age range 0-6 years old. A survey asked various "Choose all that apply" questions and I'd like to compare the age densities for each subset (Did choose option X, Did Not choose option X), to look for age ranges where significant differences appear to have occurred. The reported ages are actually discrete, since people round to the nearest month below about two-years-old and to the nearest half year after that. (I'm using R's `density` to come up with the monthly density from 0-6, to try to smooth through the many zeros that would occur for months at higher ages.) The overall age density is: ![enter image description here](http://i.stack.imgur.com/dXHvl.png) My first attempt was to take the ages from each partition (Checked, Not- checked), calculate a density, and do a Chi-squared (or K-L divergence or Hellinger distance) test between the two subset densities to determine statistical significance. If the differences between the two densities seem statistically and (non-subject-matter-expert) actually significant, I plotted the two densities on a single graph. For example: ![enter image description here](http://i.stack.imgur.com/Xft3T.png) But then I thought that there is an overall age density for the total population and perhaps I should "normalize" the two subset densities by dividing each by the overall density. Or maybe I should subtract the overall density from each. Those would show more of a "deviation from the whole" kind of graph. (In the example, the TRUE's much more closely reflect the overall population since they make up the vast majority: 314 out of 350.) Do any of these three ideas (raw, divide by overall density, subtract overall density) make any (statistical) sense?


title: Finding the variance of the estimator for the maximum likelihood for the Poisson distribution

text: If $K_1, \dots, K_n$ are i.i.d. Poisson distributions with parameter $\beta$ I have worked out that the maximum likelihood estimate is $$\hat\beta (k_1, \dots, k_n) = \frac{1}{n} \sum_{i=1}^n k_i$$ for data $k_1, \dots, k_n$. Therefore we can define the corresponding estimator $$T = \frac{1}{n} \sum_{i=1}^n K_i .$$ My question is how would you work out the variance of this estimator? In particular, as each $K_i$ follows a Poisson distribution with parameter $\beta$ I know, from the properties of the Poisson, that the distribution $\sum_{i=1}^n K_i$ will follow a Poisson distribution with parameter $n \beta$, but what is the distribution of $T$?


title: Bivariate conditional probability confidence

text: I have two random variables following normal distributions, $X\sim N(\mu_{x},\sigma_{x})$ and $Y\sim N(\mu_{y},\sigma_{y})$ I know the covariance is $\operatorname{Cov}(X,Y) = c$. And I believe (correct me if I'm wrong) that I can estimate the expected value of $Y$ given $X=\tau$ as follows: $E(Y|X=\tau) = \rho \frac{\sigma_{y}}{\sigma_{x}}(\tau -\mu_{x})+\mu_{y}$ where $\rho = c/(\sigma_{x}\sigma_{y})$. I believe this is a decent estimate. But now I want some way of estimating our "certainty" that our calculated value is correct. Obviously for $|\rho| \sim 1$ we can be very confident. For $\rho \sim 0.5$ a little bit confident. And for $\rho \sim 0$ I suppose our confidence isn't changed at all. i.e. since no correlation exists we have not changed any ideas about how confident we are about our expected $Y$ value. How can I express this numerically?


title: Relation between true positive, false positive, false negative and true negative

text: I have values for `True Positive (TP)` and `False Negative (FN)` as follows:               TP = 0.25     FN = 0.75      From those values, can we calculate `False Positive (FP)` and `True Negative (TN)`?


title: Cox Proportional Hazards - Dealing with unknown values for factors for some individuals in a population

text: I'm working on a model with multiple covariates in R. Unfortunately, for one of the factors in the model, the data is spotty and I have a large portion of unknown values. Rouhgly around 30% are NA's. Is it recommended to simply remove these individuals from the model (I believe this is what R does by default, "10,126 observations deleted due to missingness") or to include an Unknown category? An alternative question, is there a method to recover some of the information they yield for other factors? Thanks for any advice. If anyone knows of a good source that discusses this issue, I'd greatly appreciate it.


title: What test of independence of sample proportions is appropriate when n is very large and p is very small?

text: I have the (not necessarily veridical) impression that typical textbook independence tests of proportions are intended for use the proportions involved are not super small. In my case n may be 10^8 but the count of successes may be 10^4. N per group may also differ, which of course means the sample variances for sample proportions differ. How can I test for the independence of two sample proportions when sample size is large and the sample proportions are very small?


title: Relating conditional probability and conditional entropy

text: I am teaching undergrad computer security and look for a semi-formal way to make students understand the notion of a one-time pad (OTP). To this end, consider the following exercise: > Suppose you want to encrypt a message $M \in \{0, 1, 2\}$ using a shared > random key $K \in \{0, 1, 2\}$. Suppose you do this by representing $K$ and > $M$ using two bits (00, 01, or 10), and then XORing the two representations. > Does this scheme have the same security guarantees of the one-time pad? > Explain. The idea here is to get student to write down and analyze the following table:                    K  |  M   | E(K,M)       ------+------+--------      1  00  |  00  |  00      2  01  |  00  |  01      3  10  |  00  |  10         ----+------+----      4  00  |  01  |  01      5  01  |  01  |  00      6  10  |  01  |  11         ----+------+----      7  00  |  10  |  10      8  01  |  10  |  11      9  10  |  10  |  00      Here, it becomes obvious that some outcomes exclude certain inputs. For example, given $E(K,M) = 11$ an attacker knows that the sent message $M$ is not 0. The way I would like to teach it to students uses conditional probability: $$M\mid C = 00 \sim Unif\{0,1,2\}$$ Whereas: $$M\mid C = 11 \sim Unif\{1,2\}$$ Since the two distributions no longer exhibit the same image, the scheme exhibits a weakness. In other words, we have established that the scheme violates the independence: $$P(M \mid C) = P(M) \quad \forall C$$ My question is now, how can I relate this factoid to the OTP property that possession of the ciphertext does not reveal anything about the plaintext, in terms of conditional entropy: $$H(M \mid C) = H(M)$$ I am looking for an easy mathematical link between the conditional entropy and conditional probability, since only the latter represents a concept that students are familiar with.


title: Consistency of sample quantile

text: I'm new to statistics and recently I learnt about non-parametric estimation. On the estimation of empirical quantiles, many notes talk about the asymptotic behavior of sample quantiles. The proof is based on the consistency of the empirical quantile estimator, then applying Slutsky's theorem in combination with the delta method (assuming the cdf is $C^{1}$ at the considered point). However, it seems to give a very old and trivial result stating that the **empirical quantile converges in probability to the true value** (with $C^1$ cdf hypothesis). I've asked my professor and he just give the suggestion using Hoeffding's inequality (he doesn't want to spend his time explaining such an easy problem). So, could someone please provide me the proof of this, or a reference covering this issue (e.g. a note on the internet, or link to a pdf file.)


title: Time series data to measure the effect of an intervention on inpatient utilization

text: I have a general knowledge of regression, but I don't know what is the best approach to analyze these data. The dataset represents quarterly rates of inpatient use for three provider groups `prov1`, `prov2` and `prov3`. The rates were measured on quarterly basis from `qrt-4 2005` until `qrt-3 2012`. On the 3rd qrt of 2008 an intervention was implemented by provider groups `1` and `2`, provider `3` represents all other providers of a network that did not participate of the intervention; so this is our control group. The purpose of this project is to measure the effect of the intervention between the before and after periods and between the participating and non-participating providers. Here is an example of the data:               Qrt-Yr  Period  Prov1   prov2   Prov3     Q4-2005 Before  103.4   201.1   189.9     Q1-2006 Before  175.4   224.2   184.5     Q2-2006 Before  140.2   237.3   187.7     Q3-2006 Before  196.8   255.6   183.5     Q4-2006 Before  222.1   204.7   187.7     Q1-2007 Before  172.5   212.5   189.7     Q2-2007 Before  172.7   207.9   185.9     Q3-2007 Before  128.8   214.7   192.9     Q4-2007 Before  164.9   197.7   205.5     Q1-2008 Before  172.1   225.9   210.3     Q2-2008 Before  103.9   209.3   203.1     Q3-2008 Before  151.1   223.7   211.4     Q4-2008 After   200.8   234.3   189.6     Q1-2009 After   158.1   185.6   206.8     Q2-2009 After   103.8   218.2   187.9     Q3-2009 After   132.7   204.6   202.5     Q4-2009 After   178.7   218.7   191.3     Q1-2010 After   181.0   223.6   178.4     Q2-2010 After   119.9   203.2   191.3     Q3-2010 After   144.6   228.7   193.8     Q4-2010 After   87.2    161.0   173.1     Q1-2011 After   102.4   178.7   175.8     Q2-2011 After   151.0   173.4   169.9     Q3-2011 After   137.1   233.9   187.5     Q4-2011 After   83.0    226.6   164.3     Q1-2012 After   147.6   193.7   172.0     Q2-2012 After   103.8   163.5   171.7     Q3-2012 After   159.1   180.8   172.5     Q4-2012 After   133.1   155.1   160.7


title: Skew in both directions and dealing with outliers

text: I have a lot of wonderfully messy data (got to love the social sciences), and realized that I was not fully prepared to bear its wrath. For the record, after reading some articles regarding ANOVA, I am extremely paranoid about violating its assumptions. So I am trying to do my best in cleaning the data. I have also, in instances where assumptions have been violated, also run a nonparametric counterpart to assess whether the results are similar. I'm not entirely sure whether this is an acceptable practice? Or, if so, how this can be reported in my results. A bit of background: my design is set up for a 2 x 5 mixed factorial (the between group is age - young/old - and the within groups factor is different forms of input methods). Questions:   1. In looking at my data (especially preference data) for the conditions, most of the data is skewed. The problem is that most of the distributions are positively skewed, while one or two are negatively skewed. Not quite sure what to do with this? Everything I have read has suggested that ANOVA is not a viable option under this circumstance. And some have also suggested that nonparametric tests are also not very powerful in this case. What to do?   2. I have had to transform my data to help with normality and variance issues, but still have outliers in some cases. (I have looked at the outliers and reasoned through their existence. Some people just performed very well/poorly or possibly were harsher raters than the average joe). I was just taking the extreme scores that fell above/below my threshold (defined below) and recoded them as one more (or less) the next most extreme score. My questions are: is it OK to recode outliers after the transformation? If so, what's the best way to do that, say, after a square root transformation (obviously making it one more or one less than the next most extreme score is not the way to go as with untransformed data)?   3. Originally, I started identifying outliers (in addition to normality plots and boxplots) as +/- 2.24 SD from the mean. After doing some reading, I thing that discovering outliers using the median absolute deviation might be best for my distributions. Just to be clear, though, after finding the median of the deviations of the historical scores from the median for that condition, do you just multiply that by 3 and consider all scores greater than/less than that value an outlier? Sorry, these are probably basic questions. I appreciate any feedback/advice! Thank you! **EDIT:** My dependent variables are different. Most of them are continuous (performance in WPM, and also likert-type response scales). Another type of variable is error rates (percentages). **EDIT #2:** My design is 2 (between groups - age) x 5 (within groups - input method type). There are 25 people in each group (young/old) so there are 25 people in each cell in my design. It was counterbalanced for order effects. To be more specific, my dependent variables are: 1. Adjusted Words per Minute (this is an average for each person over 15 trials) for each of 5 input methods. Along these lines, I also captured 3 different error rates (also averaged over 15 trials for each person). 2. Subjective ratings for each of the input methods: satisfaction (using the System Usability Scale, typical in our field, which is continuous, a composite score derived from 10 questions), workload (using the NASA-Task Load Index, again, continuous), pre- and post- test preference (0-50 point scale), post-test perceptions of accuracy (same 0 - 50 point scale). Also the Attitudes Toward Computers Questionnaire to assess attitudes toward technology (composites). 3. Individual measures (which I plan to correlate to some of the preference/performance measures). Individual baseline measures are all continuous, either motor speed (finger tapping test), speech rate, cognitive tests (Digit Symbol Substitution, etc.). Also some anthropometric measurements (finger circumference, hand breadth, etc.). Many of those are just going to be descriptives about my sample.


title: Unusual Panel Estimation

text: Good day! I'm using R and I have a panel dataset of 72 countries x 19 years x 6 parameters. The current mini-goal is to estimate following equation: ![enter image description here](http://i.stack.imgur.com/0Yit8.png) The problem is that _beta_ doesn't have an _i_ index, so the model cannot be estimated directly through pvcm() function of _plm_ package. Moreover, i need to estimate this equation 5 times (GDP ~ one of the parameters every time ), and in reverse (every parameter ~ GDP), so some automatization of the process would be very helpful. Thanks in advance! On second thought #1 Frisch–Waugh–Lovell theorem can somehow be used, but i have currently no idea in which way it should be done. The problems are accounting the alhpa and equation system specification if we estimate with p>1.


title: Covariance between factors in CFA (Amos)

text: I ran a Confirmatory Factor Analysis (CFA) in Amos and obtained a model with high covariance between the factors. James described this covariance in his video as alarming. What can I do to tackle this problem? All my factor loadings are high and it seems to be a good model, except that I have high covariance between the factors. Additionally, the model fit seems to be poor. All indices indicate that. Somehow, I was really surprised, because I used the TAM model and the standard items. Therefore, content validity should be given and I expected that construct validity would be given as well.


title: how is the "best" split at node determined in matlab RF package?

text: Does someone know how the "best" split at a node is determined in the MATLAB RF package? I have more than two classes, and just wondered how it is working. Is the package using the Gini impurity? (Well I would not know exactly how that works, but at least I would know then what to google for.)


title: How to test the robustness and the performance of a novel classification algorithm

text: Suppose you have a new algorithm that you want to publish. Are there any best practices and methodologies you usually consider in order to test the robustness and performance of a new method? The case is a novel classification algorithm of matrices composed by observations x features (e.g. a gene expression matrix). Thanks


title: using Root Mean Squared Error (RMSE) to compare models with different sample size

text: I'm using k-fold cross-validation to compare different models. I splitted my dataset in 6 chunks and used 4 random chunks as training set and the remaining 2 as a test set. Now I fitted n-different models to the training set and calculated the RMSE on both the training and the test sets. From what I understand, the model having the lower RMSE in the test set should be the preferable one. For the sake of clarity for I mean: RMSE = sqrt( (fitted-observed)^2/ n.observations ) The models differ one another for some indipendent variables which have differenent amounts of NA values (in particular since some variables represent the cumulative effect of others I have that the number of NAs increases the more variables I cumulate). So I find myself comparing a first model with say n NAs with a second one having 10n NAs. In this way I'm comparing models that are fitted to a different number of observations. 1) Is this an issue when comparing RMSE calculated on the test set? I know for example that, if I was comparing models on the training set, the AIC would not be meaningful in this case, less sure for the R-squared... 2) since I run each model 10 times on 10 training sets and tested on 10 test sets (see beginning for explanations), for a given model I have average RMSE and its standard error on both training and test sets. How should I interpret differences between the training and test RMSE? Any suggestion appreciated!


title: What is the purpose of characteristic functions?

text: I'm hoping that someone can explain, in layman's terms, what a characteristic function is and how it is used in practice. I've read that it is the Fourier transform of the pdf, so I guess I know _what_ it is, but I still don't understand its purpose. If someone could provide an intuitive description of its purpose and perhaps an example of how it is typically used, that would be fantastic! Just one last note: I have seen the Wikipedia page, but am apparently too dense to understand what is going on. What I'm looking for is an explanation that someone not immersed in the wonders of probability theory, say a computer scientist, could understand.


title: Treating 'Don't know/Refused' levels of categorical variables

text: I am modeling Diabetes Prediction using Logistic Regression. The dataset used is the Behavioral Risk Factor Surveillance System (BRFSS) of the Center for Disease Control (CDC). One of the independent variables is High Blood Pressure. It is categorical with the following levels 'Yes', 'No', 'Don't know/Refused'. Should I remove those rows with 'Don't know/Refused' while building the model? What difference does it make to keep or remove those rows from the model?


title: Are two identical time series cointegrated?

text: I wonder if two identical time series are cointegrated. Can anyone shed some light on this? Thanks


title: Displaying Log Periodogram

text: I want to perform spectral analysis and was wondering if this command in R:               out=mvspec(data,detrend=TRUE,log="yes")       displays the detrended log-periodogram. Also would the cepstrum be the frequency with the highest amplitude?


title: Using known correlations to select a sample that maximizes model prediction

text: I want to model the outcome of a binary dependent variable D over a large universe [many records, 100s of variables]. Given a fixed sample size, I want to maximize the accuracy of the model over the universe. I already have some information about the nature of D that I want to use. When binary variable V = 1 then D = 1 99% of the time. When V = 0 then D=1 50% of the time. Let's say half of our universe has V = 1 and the other half has V - 0. What is the optimal way to draw a random sample so that our model of D has highest accuracy? The point here is that I don't need to sample as many observations where V = 1, since that is much more easy to predict. Any insight would be most helpful!


title: Closed form formula for distribution function including skewness and kurtosis?

text: Is there such a formula? Given a set of data for which the mean, variance, skewness and kurtosis is known, or can be measured, is there a single formula which can be used to calculate the probability density of a value assumed to come from the aforementioned data?


title: Difference-in-differences estimator

text: Consider the following model:   * $wage_i=\beta_0+\beta_1after_i+\beta_2female_i+\beta_3after_ifemale_i+u_i$ where   * $after_i=1$ if date after gender-wage discrimination policy; $0$ if date before gender-wage discrimination policy.   * $female_i=1$ if female; $0$ if male. The 'treatment group' is females. This model measures the effect of a gender-wage discrimination policy on the average wage of women relative to men. The affect of the policy is captured by $\beta_3$ which can be estimated as follows:   * $\hat{\beta_3}=(E[wage_i|after_i=female_i=1]-E[wage_i|after_i=1, female_i=0])-(E[wage_i|after_i=0, female_i=1]-E[wage_i|after_i=female_i=0])$ which is known 'difference-in-differences' estimator. Questions:   1. If $\hat{\beta_3}>0$, then does this means that the policy caused women's earning to increase relative to mens?   2. For $\hat{\beta_3}$ to be interpreted as the causal effect of the policy on wages, does $E(u_i|female_i, after_i)=0$ have to be assumed? Any guidance would be very much appreciated. Thank-you.


title: How to discriminate between non-nested models?

text: ![Question](http://i.stack.imgur.com/9RaBO.png) So first this question asks "how I would discriminate between model(1) and model(2)". It appears that both models are non-nested so I would come up with a hybrid model consisting of Xt,Zt and Qt regressors. I know I need to test for Zt=0 and Qt=0 but I'm not sure if I need to test for Xt=0. My logic tells me no as the point is to compare model(1) and model(2), since they both contain Xt as regressor this won't tell me if model(1) is better than model(2) and vice- versa. But someone please correct me if I'm wrong. For the second part where: Given the available information, would it be desirable for you to come up with a third and dominant model? First a third model would be desirable only if model(1) and model(2) are not "good" or we would want a better model. But from the way the question is formulated it cannot be the second one, so I believe I should test each individual models to see if they are "good". For this part the only thing I can think about is doing the following for each equation: t-test(Individual Significance Test) F-test(Joint Significance Test) Apart from these two tests there aren't really any other information that can be used to tests the models. So if model(1) and model(2) fail any of these tests then a third and dominant model would be desirable. Is my understanding correct?