title: R scatterplot matrix with nonparametric density

text: I normally use MATLAB, or JMP but right now am working with R. I have ~150 dimensional data with a few hundred thousand rows. Some of the columns are non-informative, they only have one value. This makes some of the descendants of "pairs" to fail. I wish that I could say "pairs(mydata)" but that doesn't work well. I know a few (~18) higher value variables that I think are more informative so I can, if I like, enter "pairs(mydata[,indices])" and get an 18x18 scatterplot matrix plot. I don't get the nonparametric density on top of them. Can you tell me how to get this? Now JMP has a switch that allows adding what they call a nonparametric density to scatterplots. (link, link,) I want to do that for my data. Here are things that I have looked out but they have not been helpful:   * (pairs) http://www.statmethods.net/graphs/scatterplot.html   * (ggpairs) http://www.r-bloggers.com/five-ways-to-visualize-your-pairwise-comparisons/   * http://learnr.wordpress.com/2009/07/15/ggplot2-version-of-figures-in-lattice-multivariate-data-visualization-with-r-part-5/ I really want something like this, with my data, as each of the scatterplot panes in the scatterplot matrix: ![enter image description here](http://i.stack.imgur.com/qTsca.png)


title: In testing, do we need to make the area of an acceptance region as small as possible?

text: One criterion to access a confidence interval is that the smaller its area, the better. For hypothesis testing, I was wondering if it is also that the smaller the area of the acceptance region, the better? When we pursue maximum power while restricting type I error to be no more than the significance level, is the area of the resulted acceptance region usually or necessarily small(est)? Thanks and regards!


title: How to test for one way correlation between time series?

text: I want to test existence of correlation between two time series but I suspect that there would exist only one way correlation. By one way, I mean that changes in Series A should bring about a change in Series B, but the converse may not be true. Do any tests cater to this situation?


title: Should I use aggregate or disaggregate data?

text: I have this dataset of flows:               Source entity | Dest entity | Traffic | Cost | Source location | Dest location | Direction | more independent variables (mostly nominal)      There is also a unit price ($/unit traffic) of each entity which comes from a discrete set `{p1,p2,p3}` and is ordinal/continuous. I want to model this unit price using regression analysis. Now the question i'm facing is that price is assigned to an entity (which can be source or dest in the table above) and not to a flow which represents each row above. I'm assuming that flows are independent of each other (i.i.d). Would it be wise to attribute unit price of an entity to a flow ? I know there is the option of aggregating data on entities but i'm afraid that could be disadvantageous because:   1. It is likely that information would be lost   2. Dependency is introduced among rows Also, which models could make sense ? I'm inclined towards regression because of simplicity. Appreciate any help/references here. Thanks PS. I'm already confused while dealing with this many nominal variables.


title: How to choose interactions to test when there are many predictors?

text: If I have many covariates $X_1, X_2, ... ,X_{15}$ in a linear regression model, how do I determine which two-way interactions to include? Obviously there would be too many potential models to do a best subset regression. Should I fit all the two-way interactions and do a stepwise selection?


title: Why is ANOVA interaction no longer significant after including a covariate?

text: I have initially run a 2x2x2 Anova and found some significant interaction between the independent variables. I then ran a ANCOVA where I added a covariate. I am confused with how to interpret the output from the ANCOVA as the significant interactions that I found with the ANOVA are now insignificant in the ANCOVA.


title: How do I setup a model with hierarchical structure using lmer in R?

text: I am trying to isolate the important predictors for my response variable "Y". I know that "TL" (which is an individual level predictor) affects "Y", and now I want to determine if adding the site level predictors "PC1" and "PC4" (both continuous variables) improve the fit of the model. The data looks like this:               Y    TL    PC1    PC4    Site     5.6  17    -2.26  1.89   Site A     5.9  19    -2.26  1.89   Site A     5.2  20    -2.26  1.89   Site A     4.0  17    -1.56  2.34   Site B     5.0  19    -1.56  2.34   Site B     5.6  18    -1.56  2.34   Site B      Overall I have approx. 15 samples at thirteen sites. Here are three models I setup:               model1<-lmer(Y~ TL + (1 + TL | Site), data=sauro.all) #(no site level predictors)     model2<-lmer(Y ~ TL + PC1 + (1 + TL | Site), data=sauro.all) #PC1 as site predictor)     model3<-lmer(Y ~ TL + PC4 + (1 + TL | Site), data=sauro.all) #PC4 as site level predictor      Summary of Model 1:               Linear mixed model fit by maximum likelihood      Formula: Y ~ TL + (1 + TL | Site)         Data: sauro.all         AIC   BIC logLik deviance REMLdev      400.5 419.6 -194.3    388.5   395.4     Random effects:      Groups   Name        Variance   Std.Dev. Corr        Site     (Intercept) 0.24335394 0.493309                      TL          0.00043485 0.020853 1.000       Residual             0.40595113 0.637143            Number of obs: 178, groups: Site, 13          Fixed effects:                  Estimate Std. Error t value     (Intercept) -16.84997    0.41581  -40.52     TL            0.10890    0.02012    5.41          Correlation of Fixed Effects:        (Intr)     TL -0.803      Summary of Model 2:               Linear mixed model fit by maximum likelihood      Formula: Y ~ TL + PC1 + (1 | Site)         Data: sauro.all         AIC   BIC logLik deviance REMLdev      396.3 412.2 -193.1    386.3   395.1     Random effects:      Groups   Name        Variance Std.Dev.      Site     (Intercept) 0.67023  0.81868       Residual             0.40740  0.63828      Number of obs: 178, groups: Site, 13          Fixed effects:                  Estimate Std. Error t value     (Intercept) -16.89643    0.45413  -37.21     TL            0.11022    0.01927    5.72     PC1          -0.32754    0.17083   -1.92          Correlation of Fixed Effects:         (Intr) TL         TL  -0.859            PC1  0.054 -0.054      ...and summary of model3:               Linear mixed model fit by maximum likelihood      Formula: Y ~ TL + PC4 + (1 + TL | Site)         Data: sauro.all       AIC   BIC logLik deviance REMLdev      399 421.2 -192.5      385   393.3     Random effects:      Groups   Name        Variance   Std.Dev. Corr        Site     (Intercept) 0.13411114 0.366212                      TL          0.00044998 0.021213 1.000       Residual             0.40582841 0.637047            Number of obs: 178, groups: Site, 13          Fixed effects:                  Estimate Std. Error t value     (Intercept) -16.83238    0.40527  -41.53     TL            0.11003    0.02016    5.46     PC4          -0.45466    0.22571   -2.01          Correlation of Fixed Effects:         (Intr) TL         TL  -0.845            PC4  0.005 -0.052      Model 2 has the lowest AIC score so I think this is the best model. I am having trouble with interpretation. My questions are: *How do I interpret the contribution of PC1 to model2 without p values? *This analysis is telling me how PC1 affects intercepts of (Y~TL) among sites, correct? However, slopes of this relationship could also vary with PC1. How do I check this? *What is the best way to visualize an lmer model? Thank you!


title: How to test a mediation model in panel data?

text: I have a research design with a mediation effect between 3 variables over a period of 5 years (panel data). And I would appreciate it if you can provide advice for analysis. **the independent variable** = CEO which include the eduation(level of education: PHD, MSC, etc), workexperience info (work exp in years & functions), & demographic details ( gender, age) **Dependent variable** = FIRM performance (market value, sales, net income) **Mediation variable** = diversity of partners (total amount of partners a company has, in which industry'multiple' and which function'multiple'). **Panel Data** Have collected all data over a 5-year period. (changes in all 3 variables) I have used 0,1 to for CEO work experience, function and industry I have used 0,1 for partners as well, both for the industries and functions. After consulting with some people I was suggested to use STATA to analyse the paneldata & relationships over the 5 year period. Can you please provide suggestions for the analysis. Thank you very much. C.A.


title: Is spectral leakage a random variable?

text: Spectral leakage is obviously a function of several factors, including its underlying algorithm, precision of variables, etc. Is there a body of work that explores the behavior and attempts to characterize the noise in spectral leakage? Both frequentist and Bayesian approaches?


title: Cluster or factor analysis?

text: I read a lot of forum to understand this, but I'm only more confused. I have a database with some comorbidities and I want to see if they could be divided in groups. I did a cluster analysis and a factor analysis and the results are similar. The difference is only in one comorbidity that is in another group. After I want to see if these groups have different role in functional independence. Only one article in literature used the factor analysis method. Which is the correct method I have to use?


title: LibSVM and matlab for time series prediction

text: I have a problem making time series predictions with SVM and Matlab. I tried to solve the problem by myself in several ways without success. I downloaded, compiled and installed LibSVM scripts for Matlab. But I don't know how to format my data in input. I mean, I have a time series that is an array of values, something like:               x=[1 2 3 4 5 6 7 8 9 10]      I want to use 70% of the vector as a training sequence of the SVM:               x1=[1 2 3 4 5 6 7]      Then I have to predict the last three values of the time series and I have to calculate the error between predicted values and the array `[8 9 10]`. I have to understand:   1. How to set SVM parameters. But there is a script that I found here that is the Matlab version of easy.py. So problem solved.   2. How to format my time series in order to be "acceptable" with easy.py (easy.m), svmtrain and svmpredict. On LibSVM FAQ website I found this procedure [code]:               matlab> SPECTF = csvread('SPECTF.train'); % read a csv file     matlab> labels = SPECTF(:, 1); % labels from the 1st column     matlab> features = SPECTF(:, 2:end);     matlab> features_sparse = sparse(features); % features must be in a sparse matrix     matlab> libsvmwrite('SPECTFlibsvm.train', labels, features_sparse);      But I think this is a procedure for classification. From my point of view, how do I distinguish labels and features if I have only a simple array of values? Here I found a similar procedure, but it is not clear if the script really works or not (I obtain a translated version of the input) Can anyone help me?


title: Fisher's exact test in R

text: First, I am no expert but I am analyzing some marketing data. I have information on two versions of the same site, and I have data on the number of times people filled out a form on each version of the site. I want to know if one of the site variation performs better at generating more filled out forms. Sample data:               dat2 = matrix(c(10,50,35,40), ncol=2)     dat2                                  Site 1                Site 2     Filled out form                10                    35          Did not fill out form          50                    40          > fisher.test(dat2)              Fisher's Exact Test for Count Data          data:  dat2      p-value = 0.0002381     alternative hypothesis: true odds ratio is not equal to 1      95 percent confidence interval:      0.09056509 0.54780215      sample estimates:     odds ratio       0.2311144       I'm really not sure if I set up the test properly, but I can obviously reject the null hypothesis given the low p-value. Site 2 converts better than site 1 at a statistically significant threshold. Given the problem, am I performing the correct test?


title: Meta-regression with dependent variables of different distributions

text: The textbook examples of metaregression I demonstrate how to combine coefficients when the dependent variable has the one distribution (for instance, exclusively binomial or exclusively continuous). Is there a method of combining regression coefficients for a meta-regression when the dependent variables are a mix of distributions: for instance, counts, ordinal and continuous?


title: MATLAB implementaion of kernel ridge regression

text: I am trying to implement the Kernel Ridge Regression algorithm but I am getting some quite strange results. I am afraid that I have made some silly mistakes, so I need your help to find out how to fix them. I quote my code bellow. All suggestions are welcome.               function [MSE_train_mean,MSE_train_std,MSE_test_mean,MSE_test_std,w_star] = perform_krr(X,Y)          % Set number of runs.     runs = 2;          % Set number of folds.     folds = 5;          % Set values and number of powers.     power = 0:9;     powers = size(power,2);          % Set values and number of gammas.     gamma = 2 .^ (-power);     gammas = size(gamma,2);          % Set values and number of sigmas.     sigma = 2 .^ (power);     sigmas = size(sigma,2);          % Initialize auxiliary matrices.     MSE_train = zeros(folds,sigmas,gammas);     MSE_train_optimum = zeros(runs,1);     MSE_test = zeros(folds,sigmas,gammas);     MSE_test_optimum = zeros(runs,1);     w_star_optimum = zeros(runs,round(size(X,1)*4/5));          % Perform runs.     for r=1:runs,          % Split original dataset into training and test set.     split_assignments = cross_val_split(folds,size(X,1));          % Perform folds.     for f=1:folds,          % Assign explanatory variables (x).     x_train = X((split_assignments(:,f)==0),:);     x_test = X((split_assignments(:,f)==1),:);          % Assign respond variable (y).     y_train = Y((split_assignments(:,f)==0),:);     y_test = Y((split_assignments(:,f)==1),:);          % Retrieve size of matrices.     [l_train_n,l_train_m] = size(x_train);     [l_test_n,l_test_m] = size(x_test);          % Perform sigmas.     for s=1:sigmas,          % Construct the Gaussian Kernel for given sigma.     [kernel_train] = calculate_krr_gaussiankernel(x_train,x_train,sigma(s));          % Perform gammas.     for g=1:gammas,          % Compute the Dual Weights for given gamma.     a_star = calculate_krr_dualversion(kernel_train,y_train,gamma(g));          % Apply the learned weights on training set and compute the corresponding     % MSE for the given values of sigma and gamma.     MSE_train(f,s,g) = calculate_krr_dualcost(kernel_train,y_train,a_star);          end     end     end          % Average over folds.     MSE_train_mean = mean(MSE_train,1);     MSE_train_std = std(MSE_train,0,1);          % Reshape both mean and std datasets.     MSE_train_mean = reshape(MSE_train_mean,[sigmas gammas]);     MSE_train_std = reshape(MSE_train_std,[sigmas gammas]);          % Plot mean and std of training set MSEs as a function of sigma and gamma.     figure, mesh(MSE_train_mean), title('Mean of Training Set MSEs vs Gamma and Sigma');     figure, mesh(MSE_train_std), title('Std of Training Set MSEs vs Gamma and Sigma');          % Choose the optimum values for the regularization parameters, given the     % fact that we want a minimal training set error;  we use the training set     % error as a guidance to select the regularization parameters.     [V,I] = min(MSE_train_mean(:));     [S,G] = ind2sub(size(MSE_train_mean),I);          % Retrieve the optimal parameters.     optimum_sigma = sigma(S);     optimum_gamma = gamma(G);          % Construct the training and test sets: let the first 4/5 entries be the     % training set and the last 1/5 entries be the test set.     x_train = X(1:l_train_n,1:l_train_m);     x_test = X(l_train_n+1:end,1:l_train_m);     y_train = Y(1:l_train_n,1:1);     y_test = Y(l_train_n+1:end,1:1);          % Construct the optimum Gaussian Kernel of both training and test sets.     [kernel_train_optimum] = calculate_krr_gaussiankernel(x_train,x_train,optimum_sigma);     [kernel_test_optimum] = calculate_krr_gaussiankernel(x_test,x_test,optimum_sigma);          % Compute the optimum Dual Weights of both training and test sets.     a_star_train_optimum = calculate_krr_dualversion(kernel_train_optimum,y_train,optimum_gamma);     a_star_test_optimum = calculate_krr_dualversion(kernel_test_optimum,y_test,optimum_gamma);     w_star_optimum(r,:) = a_star_train_optimum';          % Apply the optimum learned weights on both training and test sets and     % compute the corresponding MSEs.     MSE_train_optimum(r,1) = calculate_krr_dualcost(kernel_train_optimum,y_train,a_star_train_optimum);     MSE_test_optimum(r,1) = calculate_krr_dualcost(kernel_test_optimum,y_test,a_star_test_optimum);          end          % % Calculate the training set's (in-sample) residuals.     % y_hat_train = x_train * a_star_optimum;     % e_train = y_train - y_hat_train;     % % Check that they sum up to zero.     % fprintf('Training set`s (in-sample) residuals sum up to %i.\n',sum(e_train));     % % Finally, plot them.     % figure, scatter(1:1:l_train_n,e_train);          % % Calculate the test set's (out-of-sample) residuals.     % y_hat_test = x_test * a_star_optimum;     % e_test = y_test - y_hat_test;     % % Check that they sum up to zero.     % fprintf('Test set`s (out-of-sample) residuals sum up to %i.\n',sum(e_test));     % % Finally, plot them.     % figure, scatter(1:1:l_test_n,e_test);          % Compute mean and standard deviation of both training and test set MSEs.     MSE_train_mean = mean(MSE_train_optimum,1);     MSE_train_std = std(MSE_train_optimum,0,1);     MSE_test_mean = mean(MSE_test_optimum,1);     MSE_test_std = std(MSE_test_optimum,0,1);     w_star = mean(w_star_optimum);          % Print summary statistics.     fprintf('Over %i runs, %i folds, %i gammas, %i sigmas, and regarding the training set, the mean and standard deviation is %.2f and %.2f respectively.\n',runs,folds,gammas,sigmas,MSE_train_mean,MSE_train_std);     fprintf('Over %i runs, %i folds, %i gammas, %i sigmas, and regarding the test set, the mean and standard deviation is %.2f and %.2f respectively.\n',runs,folds,gammas,sigmas,MSE_test_mean,MSE_test_std);          end          function [K] = calculate_krr_gaussiankernel(Xi,Xj,S)         K = zeros(size(Xi,1),size(Xj,1));         for Ixi = 1:size(Xi,1),             for Ixj = 1:size(Xj,1),                 K(Ixi,Ixj) = exp((-norm(Xi(Ixi,:) - Xj(Ixj,:)) .^ 2) ./ (2 * (S .^ 2)));             end         end     end          function [A] = calculate_krr_dualversion(K,Y,G)         A = (K + G * size(Y,1) * eye(size(Y,1))) \ Y;     end          function [C] = calculate_krr_dualcost(K,Y,A)         C = (1 / size(Y,1)) * ((K * A - Y)') * (K * A - Y);     end


Let $X_1, X_2,\cdots,X_n$ be a random sample from the distribution with the p.d.f $$f(x)=\frac{1}{\beta -\alpha},\alpha<x<\beta $$ where $0<\alpha<\beta<\infty$. Obtain the minimum variance unbiased estimators of $\frac{\alpha+\beta}{2}$ and $\beta-\alpha$. Here I try to use Rao Blackwell Method but I am not able to solve with that. Please help


title: Measure of similarity or distance between two symmetric covariance matrices

text: are there any measures of similarity or distance between two symmetric covariance matrices (both having the same dimension)? I am thinking here of analogues to KL divergence of two probability distributions or the Euclidean distance between vectors except applied to matrices. I imagine there would be quite a few similarity measurements. Ideally I would also like to test the null-hypothesis that two covariance matrices are identical.


title: A problem using Bayesian modeling to impute a variable measured with error

text: I'm trying to fit a basic measurement error model (from Wansbeek and Meijer pg. 191), using a Bayesian latent variable model. The model appear to converge, but on the wrong answer. I've tried all kinds of variations, and can't get it to work. Why? Here's the setup: We wish to regress y ~ chi. Unfortunately, instead of chi we have x, where x = chi + v, where v is a zero-mean error term. Chi itself is determined in part by W, so we will use W and x to impute chi. Here's the R code I'm using to generate data:               n <- 1000         alpha <- 2     beta <- 5          W <- runif(n)     v <- runif(n)*5     u <- runif(n)     e <- runif(n)          chi <- W * alpha + u     x <- chi + v     y <- chi * beta + e      Here's the R code I'm using to set up the JAGS model:               jags.data = list("x","y","W","n")     jags.params = c("alpha", "beta", "s_y", "s_x", "s_chi" )#, "chi")          jagsfit <- jags(jags.data, inits=NULL, jags.params,                     n.iter=50000, model.file="jags_model_1.txt")      Here's the JAGS model itself:               #jags_model_1.txt          model{         for( i in 1:n ){             y[i]  ~ dnorm( chi[i] * beta, s_y )             x[i]  ~ dnorm( chi[i], s_x )             chi[i]~ dnorm( W[i] * alpha,    s_chi )         }              alpha ~ dnorm( 0, 20 )         beta ~ dnorm( 0, 20 )              s_y  ~ dunif(0,20)         s_x  ~ dunif(0,20)         s_chi~ dunif(0,20)     }      And here's the kind of output it gives me:               > print(jagsfit)     Inference for Bugs model at "jags_model_1.txt", fit using jags,      3 chains, each with 50000 iterations (first 25000 discarded), n.thin = 25      n.sims = 3000 iterations saved               mu.vect sd.vect     2.5%      25%      50%      75%    97.5%  Rhat n.eff     alpha       5.539   0.084    5.378    5.480    5.535    5.595    5.706 1.020   110     beta        2.511   0.035    2.437    2.487    2.511    2.535    2.576 1.030    89     s_chi       1.385   0.077    1.241    1.331    1.386    1.437    1.539 1.010   220     s_x         0.326   0.015    0.297    0.316    0.326    0.336    0.357 1.011   200     s_y        15.711   3.215    8.175   13.651   16.379   18.381   19.859 1.005   540     deviance 4066.079 242.908 3769.635 3887.593 4003.225 4188.925 4703.454 1.004   580          For each parameter, n.eff is a crude measure of effective sample size,     and Rhat is the potential scale reduction factor (at convergence, Rhat=1).          DIC info (using the rule, pD = var(deviance)/2)     pD = 29419.3 and DIC = 33485.4     DIC is an estimate of expected predictive error (lower deviance is better).      The problem is that the parameter estimates are nowhere near the original values. I've run this simulation several times with different settings (priors, and jags inits) and it's never come close to the true parameters. What's wrong here?


title: How should I classify stores based on the demographics of their customers?

text: I've got a dataset of demographic details of store customers and which store they (most frequently) visit. **I would like to categorize the stores based on their customers.** To clarify: The issue here is to create clusters of shops, on the basis of the characteristics of the customers who have attended them. In other words, the aim is to create clusters of shops having a similar clientele. I have around 7,000 customer records, distributed (unevenly) across about 50 stores. Most of the customer data is categorical, but there are a couple of continuous variables. How should I go about categorizing the different stores?


title: Why is the k-means algorithm minimizing the within cluster variance?

text: I have read that the k-means algorithm tries to minimize the within cluster sum of squares (or variance). With some brainstorming, a question popped up. Why is it that k-means or any other clustering algorithm that has **within cluster variance** as its objective to minimize, chose this as the objective function to minimize? What is it about **within cluster variance** that helps you decide that this is what you want to focus while clustering? - And especially clustering? Let me put the question in another way (This question can be a sub-question or another way of putting the same question). Why would you say that minimizing **within cluster variance** is the right way of clustering (referring to the algorithms that minimize it)? Can there be other objective functions that can be minimized (or maximized or anything) for clustering?


title: Correlated factors in a global model

text: I have developed three a-priori candidate models that represent factors at three separate spatial scales. A factor in one model is correlated with a factor in another model. (there are debates in the literature about which factor is most influencing the independent variable, and I do not want to draw conclusions about one factor without giving consideration to another.) But I also learned that I need to create global model that uses all factors from all models. So, is it valid for me to create a global model that has correlated variables?


title: Are matrix decomposition based Kalman filter algorithms faster or more robust?

text: I have been using linear Kalman Filters for several different applications. I wrote the implementation from scratch and it follows Welch & Bishop verbatim in the simplest way. I have also heard of implementations that employ matrix decompositions, such as SVD. I do a simple matrix inversion in the computation of the K matrix. My question is what is the advantage of using a less naive approach? Is it for computation speed, robustness or both (or something else)? I haven't had any issues with my implementation, although my data is reasonably well behaved, at least for the moment.


title: Visualizing mixed model results

text: One of the problems I've always had with mixed models is figuring out data visualizations - of the kind that could end up on a paper or poster - once one has the results. Right now, I'm working on a Poisson mixed effects model with a formula that looks something like the following: ` a <- glmer(counts ~ X + Y + Time + (Y + Time | Site) + offset(log(people))` With something fitted in glm() one could easily use the predict() to get predictions for a new data set, and build something off of that. But with output like this - how would you construct something like a plot of the rate over time with the shifts from X (and likely with a set value of Y)? I think one could predict the fit well enough just from the Fixed effects estimates, but what about the 95% CI? Is there anything else someone can think of that help visualize results? The results of the model are below:               Random effects:      Groups     Name        Variance   Std.Dev.  Corr                Site       (Intercept) 5.3678e-01 0.7326513                                time        2.4173e-05 0.0049167  0.250                         Y           4.9378e-05 0.0070270 -0.911  0.172           Fixed effects:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -8.1679391  0.1479849  -55.19  < 2e-16     X            0.4130639  0.1013899    4.07 4.62e-05     time         0.0009053  0.0012980    0.70    0.486         Y            0.0187977  0.0023531    7.99 1.37e-15          Correlation of Fixed Effects:              (Intr) Y    time     X      -0.178                   time    0.387 -0.305            Y      -0.589  0.009  0.085


title: Covariance matrix of distance matrix from an estimated squared distance matrix

text: I would like to perform fuzzy clustering on a data set $\boldsymbol{X}$ that contains missing elements. The fuzzy clustering algorithm requires the computation of covariance matrix $\boldsymbol{C}$ of the distance matrix $d_i = \boldsymbol{x}_i - v_j$ where $\boldsymbol{x}_i$ is a row vector from $\boldsymbol{X}$ and $v_j$ is the center of the jth cluster, also a row vector. Since $\boldsymbol{X}$ consists of incomplete samples, i.e., $\boldsymbol{x}_i$ might have missing element(s), $d_i$ cannot be computed. However, I can estimate the squared distance $d_i^2$ for the incomplete rows, and certainly for the complete rows. Because I need to compute the covariance matrix of the distance, is it possible to get the (estimated) covariance matrix of the distance by taking the square root the covariance matrix of the squared distance?


title: First iteration in MCMC coda chain is different from initial values

text: I have my jags output object. In order to understand how MCMC coda chains work, I tried to see if first iteration in each MCMC chain is equal to the initial values supplied. And it is different! The initial value is not there! Is it an error? Note that I specified burnin = 0 for this purpose. How I ran jags:               inits = function () { list(         alpha = rnorm(no_crit, 0, 10000),         beta = rnorm(no_crit, 0, 10000)         ,eps_tau = 7.9         ,gamma_tau = 3.1         ,delta_tau = 213      )      }          params = c("alpha", "beta", "eps_tau", "gamma_tau", "delta_tau")          ni <- 5000     nt <- 8     nb <- 0     nc <- 3          out <- R2jags::jags(win.data, inits, params, "model.txt",         nc, ni, nb, nt,           working.directory = paste(getwd(), "/tmp_bugs/", sep = "")     )      After the jags computation finished, I dumped the first iteration from each MCMC coda chain:               > mm = as.mcmc(out)     > mm[1, c("delta_tau", "eps_tau")]     > mm[1, c("delta_tau", "eps_tau")]     [[1]]        delta_tau      eps_tau      4426.7716020    0.4825011           [[2]]        delta_tau      eps_tau      4811.3174529    0.5240721           [[3]]        delta_tau      eps_tau      4406.2672016    0.5351576       As you can see, the first iteration in all these chains is different from what I supplied as initial values (eps_tau = 7.9, delta_tau = 213).


title: Interpreting R summary output

text: In an assignment, I am provided with the following output from R:               > summary(myfit)     Call:     lm(formula = Blood_pre ~ Heart_wt)          Residuals:         Min      1Q  Median  3Q         Max     -41.655 -29.796  -8.198  27.425  48.161          Coefficients:                 Estimate Std. Error t value Pr(>|t|)     (Intercept)    82.413    30.983   2.660  0.02219 *     Heart_wt        5.723     1.779   3.218  0.00819 **     ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 33.72 on 11 degrees of freedom     Multiple R-squared: 0.4849,     Adjusted R-squared: 0.4381     F-statistic: 10.35 on 1 and 11 DF, p-value: 0.00819      ![graph](http://i.stack.imgur.com/WXZ2S.png) (I am also given that $n = 13$) From that above data, I am suppose to suppose to find the equation of the least squares line, of which I know the formula. In that formula, however, there are the variables $SS_{xy}$, $SS_{xx}$, $\bar{x}$, and $\bar{y}$. Without the original dataset, I have no idea how to calculate those variables. Could someone help point me in right direction?


title: Matlab tool for nonlinear mixed effects

text: I'm wondering if someone is familiar with the NLMEFITSA algorithm in matlab. This algorithm gives me as result the fixed effects parameters (beta) for a mixed effects model as well as their covariance matrix (psi), but it also gives me, in the "stats" struct, one covariance matrix called "covb" I know this one has to do with the standard errors and it is important for calculating confidence intervals, but to be honest I don't know what is exactly the difference between this "covb" and the "psi " covariance matrix. And how could I use "covb" when simulating new data? Thanks if anyone can explain this to me.


title: Recommended books on experiment design?

text: What are the panel's recommendations for books on design of experiments? Ideally, books should be still in print or available electronically, although that may not always be feasible. If you feel moved to add a few words on what's so good about the book that would be great too. Also, aim for one book per answer so that voting can help sort the suggestions. (Community Wiki, please edit the question if you can make it better!)


title: Shifting bootstrap confidence interval to be centered around original parameter

text: I've been doing a bit of research into bootstrapping as I've been told one method of performing it, and this seems to differ from what I can find in other sources. I have a sample, and want to estimate the mean or median. I generate 1000 resamples, without replacement, calculating the mean/median for each. By taking the 26th largest and 26th smallest mean/median I calculate, I get a confidence interval. The sources I've read seem to leave it at that. However, the method I was told went one step further - it calculates the width of this interval, then generates a new interval that is centered around the mean/median of my original sample. So, if my original mean was 20, and bootstrapping gives me an interval of [17,27], it shifts this to give a final confidence interval of [15,25]. Does this make sense / have any statistical backing? Or is this a mistake?


title: Estimating effects on membership in a cluster

text: Suppose you want to find clusters based on a set of variables $Y$, and that you want to estimate the effects of some variables $X$ on membership in those clusters. Here is how I am doing it now. Step 1: Perform model-based clustering on the variables $Y$ (using the `mclust` package for this). Step 2: Optimize a multinomial regression model with cluster membership as the outcome variable. It seems like there must be a better way in which the models are estimated simultaneously. Anyone know a good tool in R for this and, even better, a good set of references for (a) the statistical model that the package implements, and (b) how to use the package? Thanks


title: OLS, HLM, or data aggregation?

text: I have some survey data from employees in different companies and I would like to perform a regression with company-level variables (based on employee perception) as both the dependent and independent variables. I am wondering which strategy to use - straight OLS, hierarchical linear modeling, or mean aggregation across organizations (e.g., the mean score for employees in different organizations). For example, the dependent variable could be 'perceived company performance' and the independent variables 'managerial quality' and 'performance bonus.' I also have a number of demographic variables (sex, tenure, etc.). I have already performed all of these analyses, and the results are very much similar (I included the demographic variables only in the OLS and HLM regressions). From a methodological point of view, which approach has the strongest theoretical foundation? Please let me know if I can provide more information. Thanks.


title: Markov chain question

text: Let S and C denote a Sunny day and a Cloudy day, respectively. Suppose that the weather on the given day depends on the past weather only through the weather on the two immediately preceding days according to the following conditional probabilities: The probability that a day is sunny given that the immediately preceding two days were sunny is P(S|S,S) = .8. The probability that a day is sunny given that the immediately preceding day was sunny and the day before was cloudy is P(S|C,S) = .6.The probability that a day is sunny given that the immediately preceding day was cloudy and the day before was sunny is P(S|S,C) = .4. Similarly, P(S|C,C) = .1, P=(C|S,S) =.2, P(C|C,S) = .4, P(C|S,C) = .6, P(C|C,C) = .9. Let n = 0 denote the initial day (day 0), let n = 1 denote the 1st day, and so on. We also know the following two initial probabilities: Prob(the initial day is cloudy, the first day is sunny) = .1 and Prob (the initial day is cloudy, the first day is cloudy) = .9. a.) What is the probability that the second and third day will be cloudy? b.) What is the probability that the third day will be cloudy? c.) What is the long-term probability that a day will be cloudy?


title: Algorithms/methods to create more features of a limited amount of features?

text: So, let's suppose that I have a set of 20 features - some of them are continous and some of them are binary. Is there an algorithm/method/solution to create more features ( combine/transform ) those features ( for example by using 2 of them to create another feature )? I want to do this to give a linear SVM more features to look at - is this a good thing to do or a big no-no? Would this be called increasing the dimensionality of the dataset? Can I extract more data out of those 20 numbers by creating new variables of their combinations? If this matters at all, I use a linear kernel.


title: What is the correct likelihood function for an sequential, adaptive data generation process?

text: Consider the following sequential, adaptive data generating process for $Y_1$, $Y_2$, $Y_3$. (By sequential I mean that we generate $Y_1$, $Y_2$, $Y_3$ in sequence and by adaptive I mean that $Y_3$ is generated depending on the observed values of $Y_1$ and $Y_2$.): $Y_1 = X_1\ \beta + \epsilon_1$ $Y_2 = X_2\ \beta + \epsilon_2$ $Y_3 = X_3\ \beta + \epsilon_3$ $ X_3 = \begin{cases} X_{31} & \mbox{if }Y_1 Y_2 \gt 0 \\\ X_{32} & \mbox{if }Y_1 Y_2 \le 0 \end{cases}$ where, $X_1$, $X_2$, $X_{31}$ and $X_{32}$ are all 1 x 2 vectors. $\beta$ is a 2 x 1 vector $\epsilon_i \sim N(0,\sigma^2)$ for $i$ = 1, 2, 3 Suppose we observe the following sequence: {$Y_1 = y_1,\ Y_2 = y_2,\ X_3 = X_{31},\ Y_3 = y_3$} and wish to estimate the parameters $\beta$ and $\sigma$. In order to write down the likelihood function note that we have four random variables: $Y_1$, $Y_2$, $X_3$ and $Y_3$. Therefore, the joint density of $Y_1$, $Y_2$, $X_3$ and $Y_3$ is given by: $f(Y_1, Y_2, X_3, Y_3 |-) = f(Y_1|-)\ f(Y_2|-)\ [\ f(Y_3|X_{31},-)\ P(X_3=X_{31}|-)$ $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad + \ f(Y_3|X_{32},-)\ P(X_3={X_{32}}|-)\ ]$ (Note: I am suppressing the dependency of the density on $\beta$ and $\sigma$.) Since the likelihood conditions on the observed data and our sequence is such that $y_1 y_2 >0$. Therefore, we have: $L(\beta,\ \sigma | X_1,\ X_2,\ X_{31}, y_1, y_2, y_3) = f(Y_1|-)\ f(Y_2|-)\ f(Y_3|X_{31},-)\ P(X_3=X_{31}) $ Is the above the correct likelihood function for this data generating process?


title: Feature selection given non-normal variables?

text: I'm trying to reduce noise (improve separability) among groups in a data set with 26 numerical variables and 10.000 samples. Each sample is a chemical profile, with each variable indicating the quantity of some chemical constituent. Most of the variables are far from normal. My question is: how can I test whether samples are from separate multivariate distributions? My intuition is to bootstrap the means for each variable within groups and then perform a permuted MANOVA test (or Hotelling's T-square) to simulate a null distribution. I know this works in the univariate case, as I've done so on microarray data, but I'm rather confused about whether or not I'm still violating assumptions. The idea is to select samples I'm confident (at some alpha) are from different distributions, then use them as the basis of a clustering function. My data is highly collinear, which may indicate that certain samples have a common origin, and therefore the assigning of novel samples to known clusters would be informative. Thanks a lot.


How can I use "partial pooling" to shrink cells with a small sample size towards the mean in a contingency table? I am interested in the proportions within each cell with a positive response. To clarify, I am tabulating various customer demographics survey across a range of stores. I would consider a standard contingency table with no shrinkage a "no pooling" estimate, and a mean proportion for all stores a "total pooling" estimate. I am interested in a compromise between the two, with stores with a small sample size having their variance and point estimates shrunk towards the total pooling estimate. This picture is further complicated because the proportions are typically quite low, with a total pooling estimate of 0.5%. What is the recommended way to construct such a contingency table? Should I use a hierarchical model with only a store level variables? And, if so, what sort of GLM should I use? Some synthetic data with similar characteristics to mine is below:               set.seed <- 2012          gen.prop <- function(n, x)sample(c(rep(T, round(n*x, 0)), rep(F, n - round(n*x, 0))), n)          df <- data.frame(Store=factor(c(rep("a", 1000), rep("b", 3000), rep("c", 1000),                                      rep("d", 20), rep("e", 10), rep("f", 10))),                       Likes.Cats=c(gen.prop(1000, 0.2), gen.prop(3000, 0.006), gen.prop(1000, 0.002),                                   gen.prop(20, 0.4), gen.prop(10, 0.1), gen.prop(10, 0)))          # Total pooling estimate     round(mean(df$Likes.Cats), 3)         # No pooling estimate         round(prop.table(table(df$Store,  df$Likes.Cats), 1), 3)


title: Identifying Condition in Conditional Probability P(A|B) or P(B|A)

text: Thanks everyone for the advice of my earlier question. May I once again get some guidance on the following? If C = { purchase by cash }, and U = { transaction amount > $200 } So if a customer purchases a new shirt that costs $80, and I need to find the probability of him paying by cash, will it be right for me to represent it as P(C|Ū)? Appreciate any advice. Many thanks in advance! =)


title: Model selection in logistic regression

text: I am building a logistic multiple regression with 5 potential variables candidates. I read this post and fit all 2^5 = 32 possible combinations of explanatory variables and chose the best model by AIC. However the ‘best’ model does not include any significant variable whereas some other models (with higher AIC) do have some significant variables. I am not statistician and do not understand why such situation. Thank you for any comments, clarification or guidance for another approach if needs be.


title: How to limit my input data for Jaccard item-item similarity calculation?

text: I'm trying to compute item-item similarity using Jaccard (specifically Tanimoto) on a large list of data in the format               (userid, itemid)      An item is considered as rated if i have a userid-itemid pair. I have about 800k users and 7900 items, and 3.57 million 'ratings'. I've restricted my data to users who have rated at least n items(usually 10). However, I'm wondering if I should place an upper limit on number of items rated. When users rate 1000 or more items, each user generates 999000 pairwise-combinations of items to use in my calc, assuming the calculation               n! / (n-r)!      Adding this much input data slows the calculating process down tremendously, even when the workload is distributed(using hadoop). I'm thinking that the users who rate many, many items are not my core users and might be diluting my similarity calculations. My gut tells me to limt the data to customers who have rated between 10 and 150-200 items but I'm not sure if there is a better way to statistically determine these boundaries. Here are some more details about my source data's distribution. Please feel free to enlighten me on any statistical terms that I might have butchered! The distribution of my users' itemCounts: ![alt text](http://www.neilkodner.com/images/littlesnapper/itemsRated.png)               > summary(raw)        itemsRated            Min.   :   1.000        1st Qu.:   1.000        Median :   1.000        Mean   :   4.466        3rd Qu.:   3.000        Max.   :2069.000            > sd(raw)     itemsRated        16.46169       If I limit my data to users who have rated at least 10 items:               > above10<-raw[raw$itemsRated>=10,]     > summary(above10)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.        10.00   13.00   19.00   34.04   35.00 2069.00      > sd(above10)     [1] 48.64679     > length(above10)     [1] 64764      If I further limit my data to users who have rated between 10 and 150 items:               > above10less150<-above10[above10<=150]     > summary(above10less150)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.        10.00   13.00   19.00   28.17   33.00  150.00      > sd(above10less150)     [1] 24.32098     > length(above10less150)     [1] 63080      Edit: I dont think this is an issue of outliers as much as the data is positively skewed.


title: X and Y are not correlated, but X is significant predictor of Y in multiple regression. What does it mean?

text: X and Y are not correlated (-.01); however, when I place X in a multiple regression predicting Y, alongside three (A, B, C) other (related) variables, X and two other variables (A, B) are significant predictors of Y. Note that the two other (A, B) variables are significantly correlated with Y outside of the regression. How should I interpret these findings? X predicts unique variance in Y, but since these are not correlated (Pearson), it is somehow difficult to interpret. I know of opposite cases (i.e., two variables are correlated but regression is not significant) and those are relatively simpler to understand from a theoretical and statistical perspective. Note that some of the predictors are quite correlated (e.g., .70) but not to the extent that I would expect substantial multicollinearity. Maybe I am mistaken, though. NOTE: I asked this question previously and it was closed. The rational was that this question is redundant with the question "How can a regression be significant yet all predictors be non-significant?". Perhaps I do not understand the other question, but I believe these are entirely separate questions, both mathematically and theoretically. My question is entirely independent from if "a regression is significant". Furthermore, several predictors are significant, while the other question entails variables not being significant, so I don't see the overlap. If these questions are redundant for reasons I do not understand, please insert a comment prior to closing this question. Also, I was hoping to message the moderator who closed the other question to avoid identical questions, but I couldn't find an option to do so.


title: Setting up this generative model for inference: uniform priors

text: I am trying to set up a generative model where I have two images $x$ and $y$ and it is assumed that $y$ can be generated by applying some unknown transformation to $x$ i.e. $$ y = t(x, w) + e $$ i.e. $$ e = y- t(x, w) $$ where $t$ is the applied transformation which has parameters $w$ and e is the error which is assumed to normally distributed IID noise. So, I have a likelihood term given by (for a single pixel in the image): $$ L(w) = \prod_{i=1}^{N}(\frac{\phi}{2\pi})^{\frac{1}{2}} \exp^{\frac{-1}{2} \phi e_i^2} $$ where $\phi$ is the noise precision. I also have a prior on the variance of $\phi$ which is modelled as a gamma distribution. Because they are independent, the likelihood terms can be multiplied as in the equation above. Now, I know absolutely nothing about the nature of the parameters $w$, so I am guessing I should have a uniform prior on $w$? However, I do not know much about how to use these priors mathematically and $w$ can take on any value, so I am not sure how to make a proper distribution out of this uniform prior as well. I was wondering if someone could give me some advice on this. I am interested in inferring maximum mode and also characterise the spread of $w$. At the moment, I am even struggling to understand how to model things in a Bayesian setting when you have no prior belief on the form of these parameters of interest.


title: some suggestions for analyzing two-factor within subject design?

text: I'm doing a two-factor within-subject design experiment, and the two factors are all categorical variables. factor one has four categories, while factor two has two categories, and is counterbalanced. this lead to a 2*4=8 conditions. the dependent variable is a continuous variable such as reaction times. my data is in a long format as follow               ID     Trial   factor1    factor2    reaction time      1       1       1          1           3.142      1       2       3          1           4.543      1       3       4          1           2.456      1       4       2          1           3.433      1       5       4          0           3.142      1       6       2          0           4.543      1       7       3          0           2.456      1       8       1          0           3.433      ...      2      2      2      ...      50      50      ...             my first question is, is it possible to run a regression analysis? (IV: two categorical variables DV: a continuous varible) my second question is , is it possible to run a regression analysis for the total participant (like 50 persons) as a whole as the same as a single person? Thanks in advance for helping!


title: How to generate nice summary table?

text: I want to have R display the data it gives me from the `summary()` function in a table so I can easily share this. I am currently just doing `summary()` in the console and then taking a screenshot, but I would rather have this generated as a nice table just like all of my graphs are. Any ideas?


title: Different parameter values when using stochastic gradient descent

text: I am having some issues with stochastic gradient descent. Using batch gradient descent where I consider all the training sets I have certain parameter values which I know are correct. My function is convex,globally. Now when I use stochastic gradient descent considering ten samples at a time, the algorithm converges nicely, but I get different parameter values. The parameter values that I get using stochastic gradient descent is like 1/2 of the one that I get using batch gradient descent. Any insights what could go wrong?


title: Flexible version of logistic regression

text: I'm trying to fit a logistic regression where there is a huge difference in the number of data points in either group (70 Vs 10,000). A statistician friend of mine has told me that this is a known problem with logistic regression and that for those kinds of numbers it overfits the data and basically doesn't work. When I bin the data and compare to the model, it is fairly obvious that this is definitely the case. I'm wondering if anyone is aware of a better / more flexible method for fitting this kind of binary response data? (I'm not a statistician by the way so go easy on me!)


title: How to do a repeated measures ANOVA with missing combinations in spss

text: I currently have a setup, for which I want do do an ANOVA in spss. In my model I need to include four factors (I'll call them `factor1` to `factor4`) each with two possible states (`true` and `false`). However not all combinations are possible. To be more precise, `factor4` can only be true if one of `factor1` or `factor2` or both are true. Hence in total I only have 14 possible combinations instead of the total 16 which would result from a full combination of all four factors. Is it possible to us a repeated measures ANOVA with such data, or do I have to use something else, since `factor4` cannot be varied independently? I tried to construct a repeated-measures ANOVA for this in the SPSS GUI, but I cannot click on OK, unless I give all 16 combinations, which is impossible given the setup I have. Is this somehow possible to program in SPSS, maybe by entering the syntax directly? Or maybe is there a way I can add the missing data, for example by adding the means of some other variable for the missing combinations, without changing the final outcome of the model? **EDIT** : I just also tried to add the missing combinations with empty variables and then run the model on this completed set of variables. However this also does not work, since in that case SPSS just complains, that there are no valid cases in the data (since all are missing those two combinations).


title: Why is generalized linear model (GLM) a semi-parametric model?

text: As we all know , the GLM has the structure: $G(EY)=X^{T}\beta$, in which $G(.)$ is a known link function. What confuse me is that some people say that it's a semiparametric model. But in my opinion, it's a parametric model, because there is no nonparametric part and all we don't know in the previuos structure is $\beta$. Could anybody tell me the reason why someone call it a semiparametric model. Thanks!


title: Logistic quantile regression – how to best convey the results

text: In a previous post I’ve wondered how to deal with EQ-5D scores. Recently I stumbled upon logistic quantile regression suggested by Bottai and McKeown that introduces an elegant way to deal with bounded outcomes. The formula is simple: $logit(y)=log(\frac{y-y_{min}}{y_{max}-y})$ To avoid log(0) and division by 0 you extend the range by a small value, $\epsilon$. This gives an environment that respects the boundaries of the score. The problem is that any $\beta$ will be in the logit scale and that makes doesn’t make any sense unless transformed back into the regular scale but that means that the $\beta$ will be non-linear. For graphing purposes this doesn’t matter but not with more $\beta$:s this will be very inconvenient. My question: **How do you suggest to report a logit $\beta$ without reporting the full span?** * * * ## Implementation example For testing the implementation I’ve written a simulation based on this basic function: $outcome=\beta_0+\beta_1* xtest^3+\beta_2*sex$ Where $\beta_0 = 0$, $\beta_1 = 0.5$ and $\beta_2 = 1$. Since there is a ceiling in scores I’ve set any outcome value above 4 and any below -1 to the max value. ### Simulate the data               set.seed(10)     intercept <- 0     beta1 <- 0.5     beta2 <- 1     n = 1000     xtest <- rnorm(n,1,1)     gender <- factor(rbinom(n, 1, .4), labels=c("Male", "Female"))     random_noise  <- runif(n, -1,1)          # Add a ceiling and a floor to simulate a bound score     fake_ceiling <- 4     fake_floor <- -1          # Just to give the graphs the same look     my_ylim <- c(fake_floor - abs(fake_floor)*.25,                   fake_ceiling + abs(fake_ceiling)*.25)     my_xlim <- c(-1.5, 3.5)          # Simulate the predictor     linpred <- intercept + beta1*xtest^3 + beta2*(gender == "Female") + random_noise     # Remove some extremes     linpred[linpred > fake_ceiling + abs(diff(range(linpred)))/2 |         linpred < fake_floor - abs(diff(range(linpred)))/2 ] <- NA     #limit the interval and give a ceiling and a floor effect similar to scores     linpred[linpred > fake_ceiling] <- fake_ceiling     linpred[linpred < fake_floor] <- fake_floor      To plot the above:               library(ggplot2)     # Just to give all the graphs the same look     my_ylim <- c(fake_floor - abs(fake_floor)*.25,                   fake_ceiling + abs(fake_ceiling)*.25)     my_xlim <- c(-1.5, 3.5)     qplot(y=linpred, x=xtest, col=gender, ylab="Outcome")      Gives this image: ![Scatterplot from simulation](http://i.stack.imgur.com/luZGu.png) ### The regressions In this section I create the regular linear regression, quantile regression (using the median) and logistic quantile regression. All estimates are based on bootstrapped values using the bootcov() function.               library(rms)          # Regular linear regression     fit_lm <- Glm(linpred~rcs(xtest, 5)+gender, x=T, y=T)     boot_fit_lm <- bootcov(fit_lm, B=500)     p <- Predict(boot_fit_lm, xtest=seq(-2.5, 3.5, by=.001), gender=c("Male", "Female"))     lm_plot <- plot.Predict(p,                   se=T,                   col.fill=c("#9999FF", "#BBBBFF"),                   xlim=my_xlim, ylim=my_ylim)          # Quantile regression regular     fit_rq <- Rq(formula(fit_lm), x=T, y=T)     boot_rq <- bootcov(fit_rq, B=500)     # A little disturbing warning:     # In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique          p <- Predict(boot_rq, xtest=seq(-2.5, 3.5, by=.001), gender=c("Male", "Female"))     rq_plot <- plot.Predict(p,                   se=T,                   col.fill=c("#9999FF", "#BBBBFF"),                   xlim=my_xlim, ylim=my_ylim)          # The logit transformations     logit_fn <- function(y, y_min, y_max, epsilon)         log((y-(y_min-epsilon))/(y_max+epsilon-y))               antilogit_fn <- function(antiy, y_min, y_max, epsilon)         (exp(antiy)*(y_max+epsilon)+y_min-epsilon)/             (1+exp(antiy))               epsilon <- .0001     y_min <- min(linpred, na.rm=T)     y_max <- max(linpred, na.rm=T)     logit_linpred <- logit_fn(linpred,                                y_min=y_min,                               y_max=y_max,                               epsilon=epsilon)          fit_rq_logit <- update(fit_rq, logit_linpred ~ .)     boot_rq_logit <- bootcov(fit_rq_logit, B=500)               p <- Predict(boot_rq_logit, xtest=seq(-2.5, 3.5, by=.001), gender=c("Male", "Female"))          # Change back to org. scale     transformed_p <- p     transformed_p$yhat <- antilogit_fn(p$yhat,                                         y_min=y_min,                                         y_max=y_max,                                         epsilon=epsilon)     transformed_p$lower <- antilogit_fn(p$lower,                                           y_min=y_min,                                          y_max=y_max,                                          epsilon=epsilon)     transformed_p$upper <- antilogit_fn(p$upper,                                           y_min=y_min,                                          y_max=y_max,                                          epsilon=epsilon)          logit_rq_plot <- plot.Predict(transformed_p,                   se=T,                   col.fill=c("#9999FF", "#BBBBFF"),                   xlim=my_xlim, ylim=my_ylim)      ### The plots To compare with the base function I’ve added this code:               library(lattice)     # Calculate the true lines     x <- seq(min(xtest), max(xtest), by=.1)     y <- beta1*x^3+intercept     y_female <- y + beta2     y[y > fake_ceiling] <- fake_ceiling     y[y < fake_floor] <- fake_floor     y_female[y_female > fake_ceiling] <- fake_ceiling     y_female[y_female < fake_floor] <- fake_floor          tr_df <- data.frame(x=x, y=y, y_female=y_female)     true_line_plot <- xyplot(y  + y_female ~ x,                               data=tr_df,                              type="l",                               xlim=my_xlim,                               ylim=my_ylim,                               ylab="Outcome",                               auto.key = list(                                text = c("Male"," Female"),                                columns=2))               # Just for making pretty graphs with the comparison plot     compareplot <- function(regr_plot, regr_title, true_plot){       print(regr_plot, position=c(0,0.5,1,1), more=T)       trellis.focus("toplevel")       panel.text(0.3, .8, regr_title, cex = 1.2, font = 2)       trellis.unfocus()       print(true_plot, position=c(0,0,1,.5), more=F)       trellis.focus("toplevel")       panel.text(0.3, .65, "True line", cex = 1.2, font = 2)       trellis.unfocus()     }          compareplot(lm_plot, "Linear regression", true_line_plot)     compareplot(rq_plot, "Quantile regression", true_line_plot)     compareplot(logit_rq_plot, "Logit - Quantile regression", true_line_plot)      ![Linear regression for bounded outcome](http://i.stack.imgur.com/74Uid.png) ![Quantile regression for bounded outcome](http://i.stack.imgur.com/xHRtF.png) ![Logistic quantile regression for bounded outcome](http://i.stack.imgur.com/XfLy8.png) ### The contrast output Now I've tried to get the contrast and it's almost "right" but it varies along the span as expected:               > contrast(boot_rq_logit, list(gender=levels(gender),      +                              xtest=c(-1:1)),      +          FUN=function(x)antilogit_fn(x, epsilon))        gender xtest Contrast   S.E.       Lower      Upper       Z      Pr(>|z|)        Male   -1    -2.5001505 0.33677523 -3.1602179 -1.84008320  -7.42 0.0000          Female -1    -1.3020162 0.29623080 -1.8826179 -0.72141450  -4.40 0.0000          Male    0    -1.3384751 0.09748767 -1.5295474 -1.14740279 -13.73 0.0000       *  Female  0    -0.1403408 0.09887240 -0.3341271  0.05344555  -1.42 0.1558          Male    1    -1.3308691 0.10810012 -1.5427414 -1.11899674 -12.31 0.0000       *  Female  1    -0.1327348 0.07605115 -0.2817923  0.01632277  -1.75 0.0809            Redundant contrasts are denoted by *          Confidence intervals are 0.95 individual intervals


title: Puzzling behavior of glmer()

text: I'd like your opinion on a very strange behavior that I recently encountered running glmer(). The problem is that when I make the dependent variable into a logical vector, glmer behaves weirdly. My dependent variable is Accuracy, and it is coded in terms of 1 (accurate response) and 0 (wrong response). What puzzles me is that transforming accuracy to a logical vector should work the same way for glmer, as a logical vector is coded in terms of TRUE or FALSE, having also 2 levels. However, glmer gives me different results depending on the transformation of the dependent variable I use. Have you guys encountered this before? Do you know why it happens? Below is sample code so you can replicate the problem yourselves.               #Create fake data     Subject   <- c(rep("S1",4), rep("S2",4), rep("S3",4), rep("S4",4))     Item      <- rep(c("I1","I2","I3","I4"),4)     Factor1   <- c(c(rep("e1",2),rep("e2",2)), c("e1","e2","e2","e1"),                 c(rep("e2",2),rep("e1",2)), c("e2","e1","e1","e2"))                       Accuracy  <- c(1,1,0,0,1,0,1,0,1,0,1,1,1,1,1,1)          #Create data frame and make "Accuracy" into a factor with 2 levels     data          <- data.frame(Subject,Item,Factor1, Accuracy)     data$Accuracy <- factor(data$Accuracy)  #Accuracy is a factor w/ 2 levels     #Run glmer     m1 <- glmer(Accuracy ~ Factor1 + (1+Factor1|Subject) + (1+Factor1|Item), family = "binomial", data= data)       summary(m1)     Fixed effects:                 Estimate Std. Error z value Pr(>|z|)       (Intercept)    1.946      1.069   1.820   0.0687 .     Factor1e2     -1.946      1.282  -1.518   0.1290       ---      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      That is the output of the first model. Now, look at what happens if I transform data$Accuracy into a logical vector when I run the model:               m2 <- glmer(as.logical(as.numeric(Accuracy)) ~ Factor1 + (1+Factor1|Subject) + (1+Factor1|Item), family = "binomial", data= data)       summary(m2)          Fixed effects:                  Estimate Std. Error z value Pr(>|z|)     (Intercept) 2.557e+01  1.259e+05       0        1     Factor1e2   2.223e-06  1.781e+05       0        1      As you can see, now the coefficient estimates are very different. As I said, this seems very puzzling to me and I'd like yo know if you have some thoughts on why this should be. Thanks a lot! \--Sol


title: Confusion relative to derivative of partition function

text: I have this partition function ![enter image description here](http://i.stack.imgur.com/ve5sd.png) Now if I take the derivative of log(Z(x)) wrt $\lambda_k$ the result is ![enter image description here](http://i.stack.imgur.com/6Qikg.png) I didn't get how this was derived. This is the paper


title: Model validation in Bayesian statistics from a model with latent variables

text: I am working with some two-regime autoregressive models first introduced by Hamilton in 1989. The specific models is of no great concern to my question, but some variables within my autoregressive models are latent binary variables - non-observable binary variables, that is. I have a dataset and use MCMC to find posterior densities of all the parameters of my models, including all the latent binary variables. I have a dataset of about 1000 observations, and for every observation I have a latent variable in my models, which can wither be 1 or 0. If my data suddendly changes, I can assume that a switch has been turned on and the latent variable has changed. So if the data, at a specific time, is (most likely) best described by the latent variable being 1, my sampler will produce a posterior distribution (for this specific latent variable) with high probability of 1 and low probability of 0. So this is all good, I have managed to do this. The question is, how do I validate my model? AIC/BIC is just for comparing models. What I am thinking is that I have to be able to plot the residuals somehow, like in normal regression where you just say that the residual is the difference of observed and predicted result. I understand that for all other parameters than the latent parameters/variables, i can just use the mean of each posterior distribution and treat it like a Maximum Likelihood Estimator. But I cant take the mean of my posterior distributions for all the latent variables, because they are either 0 or 1. It does not make any sense using, for example, 0.8 for one of my latent variables. This is no option. So how should I go on to validate my model? I am really stuck here. If I have done a bad job explaining my problem, I apologize and will try to explain better. * * * So I have been thinking. Can I use a predictive posterior distribution to validate my model? And if too many observations, lies above, say, the 95 percentile of the predictive posterior distribution, I throw away my model? Generally, for any model whose parameters are called theta, observations called x, and new observation called x_new, one calculate the predictive posterior distribution as follows: For i = 1, . . . , n (samples), we sample: 1\. theta[i] from p(theta|x) and 2\. x_new[i] from p(x_new|theta[i]) Then x_new[1], . . . ,x_new[n] are a sample from the predictive posterior distribution. But this method must surely only work when we have a iid sample. I am working with an autoregressive model. So I am stuck...


title: Sampling and Conditions

text: Suppose we consider all the people in the planet. We are interested in randomly selecting 10 people who have heart disease. Is it better to repeatedly sample 10 people and then choose the one selection in which all 10 people have heart disease? Or is it better to only look at the people who have heart disease and sample 10 from that population?


title: Determining the scale type

text: I have different functions, that measure characteristics like completeness or reputation. I try to find out on what scale types they are (nominal, ordinal, interval, ratio or absolute). Most of the functions map to values between $0.0$ and $1.0$. I would assume that they are on a ratio scale. A ratio scale is an interval scale where the values have numerical meaning, so addition and subtraction would work here. This makes somewhat sense. I represent these numbers as percentage. Is percentage in general a ratio scale? Then there are measures which map to quite arbitrarily values. For instance one to values between $-355$ and $2214$. Are these on an absolute scale type?


title: Entropy as measure of order in data streams

text: Is it sound/allowed to use entropy as a measure of (non-) uniformity of a data stream? E.g. I calculate the Shannon entropy with the standard formula based on various measures in the data stream. The measures are count values whose probabilities sum up to 1. I do the calculation continuously every 1 sec, so I get a new entropy value every second. Then I monitor the average and standard deviation of the entropy values over a sliding window, let's say the past 10 seconds. If the stddev increases or gets larger than certain threshold then I assume the stream got noiser/bumpier. Is it ok to use entropy in this context?


title: Generate random variables from the Gumbel distribution in R

text: I want to generate 1000 random numbers with Gumbel distribution (extreme value distribution type 1) with known parameters $m$ and $a$. Then I want to estimate the parameters of the distribution for these 1000 numbers using probability-weighted moments estimation (PWM) and maximum likelihood estimation (MLE) to compare them with $m$ and $a$. I appreciate your help


title: Probability when sampling from multiple urns

text: I have the following problem: I have 2 urns, one which contains which contains 20% red balls and another which contains 80% red balls. Each time I select (with replacement) a ball, I am blind and have an equal chance of selecting from each urn. **What is the probability of selecting a red ball?** It seems to me in the long run I should obtanin 50% red balls so it's 0.5. But this does seem misleading somehow? Many thanks!


title: Estimation of slope of local linear regression intercepts

text: A local linear regresson provides at each point $x_i $ an estimate of the intercept and slope of the values of $y$ centered at that point. Consider the curve formed by the sequence of intercepts. What is the 'best' estimate of its slope? Is it the estimate of the slope in the local linear regresson or is it derived from the intercepts themselves, such as might be obtained by the derivative of a spline interpolation of the intercepts. References on this topic would be appreciated.


title: Should the difference between control and treatment be modelled explicitly or implicitly?

text: Given the following experimental setup: Multiple samples are taken from a subject and each sample is treated multiple ways (including a control treatment). What is mainly interesting is the difference between the control and each treatment. I can think of two simple models for this data. With sample $i$, treatment $j$, treatment 0 being the control, let $Y_{ij}$ be the data, $\gamma_i$ be the baseline for sample $i$, $\delta_j$ be the difference for treatment $j$. The first model looks at both the control and difference: $$ Y_{ij}=\gamma_i+\delta_j+\epsilon_{ij} $$ $$ \delta_0=0 $$ Whilst the second model only looks at the difference. If we precalculate $d_{ij}$ beforehand $$ d_{ij}=Y_{ij}-Y_{i0} $$ then $$ d_{ij}=\delta_j+\varepsilon_{ij} $$ My question is what are the fundamental differences between these two setups? In particular, if the levels are meaningless in themselves and only the difference matters, is the first model doing too much and is perhaps underpowered?


title: What is compound symmetry in plain english?

text: I recently realized that a mixed-model with only subject as a random factor and the other factors as fixed factors is equivalent to an ANOVA when setting the correlational structure of the mixed model to compound symmetry. **Therefore I would like to know what does compound symmetry mean in the context of a mixed (i.e., split-plot) ANOVA, at best explained in plain English.** Besides compound symmetry `lme` offers other types of correlational structures, such as > `corSymm` general correlation matrix, with no additional structure. or different types of _spatial correlation_. **Therefore, I have the related question on which other types of correlational structures may be advisable to use in the context of designed experiments (with between- and within-subjects factors)?** It would be great if answers could point to some references for different correlational structures.


title: Sample size for binomial confidence interval

text: **Background** : A group at work is sampling 1,000 customers to contact and from that point, determining if the effort is worthwhile or not. I wanted to see if this (pretty much) arbitrary sample size value was at all "good enough". If we assume that the true proportion of successes in a population is 0.016 (1.6%) how large a sample size would I need to take in order to get a confidence interval margin of error "half width") of say 0.005 (0.5%)? Here is my approach in R:               install.packages("Hmisc")     library(Hmisc)               target.halfWidth<-0.005     sims<-25000 #number of draws from binomial to perform               p<-0.016 #true proportion     n<-seq(from=500, to=5000, by=100) #number of samples          #hold results     results<-matrix(numeric(0), length(n),2)          #loop through desired sample size options     for (i in 1: length(n))     {          x<-rbinom(sims,n[i],p) #draws from binomial with p and n      ci<-binconf(x,n[i],method="asymptotic",alpha=0.1) #normal theory 90% CI     half_width<-ci[,3]-ci[,1] #half width of CI          #need the number where the half width is within the target range     prob.halfWidth<-length(half_width[half_width<target.halfWidth])/sims          #store results     results[i,1]<-n[i]     results[i,2]<-prob.halfWidth     }          #plot     plot(results[,2],results[,1],type="b")     results       This simulation showed that we would need 2,200 samples to be 95% confident that a 90% CI would be at most 0.005. **Questions:** 1) Is this a proper method? 2) Are there better ways? 3) What advice can you give is there are finite samples of some sub- populations? Say that we want to know how many samples to take of populations where there are not "a lot" of customers to choose from. Maybe there is only 5,000 of a certain group, should we not be able to take less of them to make a determination compared to a group where there are 50,000 to choose from? Thanks!! **Add after MansT answer:** 1) Would it make sense, under my scenario with the simulation draws to add a step where this line:               prob.halfWidth<-length(half_width[half_width<target.halfWidth])/sims      only increments the numerator when the resulting CI also contains the true p (i.e. 0.016)? 2) Under your code, would it also be appropriate when dealing with a finite sample to add the finite population corrector FPC to your line:               halfWidth<-qnorm(0.95)*sqrt(p.est*(1-p.est)/n)      3) I am not sure the formula for a CI from the hypergeometric, but perhaps I could replace my code line               ci<-binconf(x,n[i],method="asymptotic",alpha=0.1) #normal theory 90% CI       with Sprop function in R?


title: Calculating confidence for a web app conversion rate

text: I'm checking the conversion rates on some segments of a web application as to figure out which types of customers we should invest in paid acquisition (theory is that these segments see more value in the product and, as such, are willing to pay). For example, on one segment there are 336 users of which 88 converted to paid. This gives a ~26% conversion rate for the segment but how can I have a high degree of confidence that the data is significant enough as to avoid burning precious resources chasing this segment? TIA


title: what kind of regression model I should use for integer and bounded dependent variable?

text: I am trying to model the rating of a product (that takes integer values between 0 and 10) using some other predictors. Can I use negative binomial regression? The data is over-dispersed toward higher ratings but I am not sure whether negbin works for bounded counts.


title: Calculating 10 year probability using hazard ratios?

text: I am working on some medical research. Given 3 values (BMI, age, gender) I have a list of baseline risks for the probability of an adverse event occurring in the next 10 years for an individual with a given BMI, age, gender. I also have a list of additional risk factors, given as Hazard Ratios. How can I properly compute the 10-year probability of an adverse event, by modifying the baseline risk if say, 1 or more risk factors apply to an individual. So, for example: given that an individual with age 30, gender Male, BMI 25 has a 5% risk of adverse event in next 10-years, assume some additional risk factors X and Y, which have respective HRs of 1.3 and 1.6, how can I modify the 5% risk to encompass the risk factors? I know that the HRs were computed using multivariable Cox proportional hazards for a 10 year period, with death considered as a competing hazard. Any additional statistical wisdom would be appreciated.


title: Generate QQ plot for sets of different size

text: I'm coming from this question Similarity between two sets of random values (which still stands BTW) where @whuber suggested in the comments that I use a Q-Q plot to assess the similarity of two discrete distributions. To recapitulate, I have two sets of random floats between [0,1] of _different sizes_ : $$A = \\{0.3637852, 0.2330702, 0.1683102, 0.2127219, 0.0152532, ..., N_A\\}$$ $$B = \\{0.4541056, 0.7521812, 0.0266602, 0.5099002, 0.3468181, ..., N_B\\}$$ where $N_B > N_A$. and I need to generate the Q-Q plot for thess sets. For what I've read, if the sets where of equal sizes ($N_B = N_A$) I would simply sort both sets from minor to major values and then plot one against the other in a 1:1 correspondence. My sets are different in sizes and the sources I've found [1], [2] both say that I need to sort both sets and then interpolate the _larger_ set ($B$ in my case). This is counter-intuitive to me since I would expect to have to generate more points from the _smaller_ $A$ set (through interpolation) so as to be able to plot this expanded set $A'$ against the $B$ set. So obviously I'm not understanding correctly the process. What are the steps I should follow to generate a Q-Q plot for two distributions like the ones I presented above?


title: Are there any known methods for picking out the "easy" cases for classification?

text: Suppose I have positive and negative samples, and I care a lot more about precision than recall. Furthermore, I assume that there are some cases among the positive samples which are "easier" to separate than others, but I have no direct way of identifying these. By easier, I mean it is easy to pick for a given classifier to pick them out without picking any negative samples. If we only had two input features, we would be able to visually pick these out as the positive examples which are visibly non-overlapping with the negative examples, but in this case I am working in high-dimensional spaces. Is there some algorithm that will, given a set of input features and labels, identify those cases which are "easy" to classify?


title: Choosing test for longitudinal, repeated measures and clustered data

text: I'm a researcher within cardiovascular medicine, currently working on a study in endocrinology (diabetes). I have the following challenge, which I have discussed with several statisticians who all gave me discordant advice and explanations (which is why I'll try the discussion here). The data set consists of:   * 400,000 observations   * 100,000 unique individuals (between 1 and 10 observations/individual).   * Varying time between observations.   * Grouping variable is ethnicity; 10 groups (basically world continents).   * A long range of covariates available.   * Outcome of interest is continuous (after log transform not satisfactory normally distributed). Thus, I have a large dataset with repeated measures with varying time between them. The outcome of interest is blood glucose levels (continuous variable). I would like to find out if/how ethnicity affects this outcome, after adjusting for covariates. The optimal statistical method for this is difficult to judge. In SAS, which I am using, there are several possibilities. PROC GENMOD does have a 'repeated' option which invokes the generalized estimated equations put forward by Liang and Zeger. This method has several advantages. PROC GLM can be used. PROC MIXED can also be used. The question is, which method is optimal for this data set and why?


title: Collinearity between categorical variables

text: There's a lot about collinearity with respect to continuous predictors but not so much that I can find on categorical predictors. I have data of this type illustrated below. The first factor is a genetic variable (allele count), the second factor is a disease category. Clearly the genes precede the disease and are a factor in showing symptoms that lead to a diagnosis. However, a regular analysis using type II or III sums of squares, as would be commonly done in psych with SPSS, misses the effect. A type I sums of squares analysis picks it up, when the appropriate order is entered because it is order dependent. Further, there are likely to be extra components to the disease process which are not related to the gene that are not well identified with type II or III, see **anova(lm1)** below vs lm2 or Anova. _Example data:_               set.seed(69)     iv1 <- sample(c(0,1,2), 150, replace=T)     iv2 <- round(iv1 + rnorm(150, 0, 1), 0)     iv2 <- ifelse(iv2<0, 0, iv2)     iv2 <- ifelse(iv2>2, 2, iv2)     dv  <- iv2 + rnorm(150, 0, 2)     iv2 <- factor(iv2, labels=c("a", "b", "c"))     df1 <- data.frame(dv, iv1, iv2)          library(car)     chisq.test(table(iv1, iv2))          # quick gene & disease relations     lm1 <- lm(dv~iv1*iv2, df1);    lm2 <- lm(dv~iv2*iv1, df1)     anova(lm1);                    anova(lm2)     Anova(lm1, type="II");         Anova(lm2, type="II")        1. **lm1** with type I SS to me seems the appropriate way to analyse the data given the background theory. Is my assumption correct?    2. I'm used to explicitly manipulated orthogonal designs, where these problems don't usually pop up. Is it difficult to convince reviewers that this is the best process (assuming point 1 is correct) in the context of an SPSS centric field?    3. And what to report in the stats section? Any extra analysis, or comments that should go in?


title: Estimating a sub-population characteristic based on independent samples without replacement

text: **Note** : This question is a follow-up of Estimating population size of a subgroup based on independent samples without replacement Let a bag have 1000 balls of arbitrary colors and unknowns sizes ($r$). Suppose we also known the total volume occupied by the balls ($t_v$). We want to estimate the total volume occupied by red balls ($t_{v}^{red}$). We attempt to do this by taking $x$ independent samples (without replacement within each sample) of size $s$ of the original population. Let $s$ be 10% of the original population as an example. Proceed to count the number of red balls in each sample, e.g., for $x = 10$:               (9, 10, 10, 11, 13, 8, 5, 15, 12, 9)      To estimate the number of red balls in the original population, one can _naïvely_ calculate the average of the above list and divide by the sample. In this case:               (9 + 10 + 10 + 11 + 13 + 8 + 5 + 15 + 12 + 9) / 10 / 0.1 = 102      What about $t_{v}^{red}$? Each sample will give us a set of sizes for the observed red balls. So:               sample1 = [3, 2, 1, 10, 2, 8, 1, 4, 1, 1] => 33 => V = 5098     sample2 = [3, 6, 1, 3, 4, 1, 1, 1, 8, 1] => 29 => V = 2673     sample3 = [1, 4, 3, 10, 1, 1, 5, 1, 1, 2] => 29 => V = 3861     sample4 = [4, 1, 1, 1, 2, 6, 2, 6, 1, 1] => 25 => V = 1624     sample5 = [1, 6, 1, 2, 4, 3, 7, 10, 10, 2] => 46 => V = 8381     sample6 = [6, 6, 1, 2, 1, 10, 3, 2, 3, 1] => 35 => V = 4728     sample7 = [1, 4, 2, 1, 4, 1, 2, 4, 10, 1] => 30 => V = 3807     sample8 = [1, 8, 3, 3, 4, 3, 1, 8, 1, 5] => 37 => V = 4074     sample9 = [7, 9, 2, 4, 4, 6, 6, 1, 8, 1] => 48 => V = 6766     sample10 = [3, 8, 1, 2, 10, 7, 3, 7, 3, 1] => 45 => V = 7191      As you can see, the estimated volume on each sample can vary considerably (from 1624 to 8381). And here lies my problem. My questions are thus:   1. How to estimate the total volume occupied by balls of a certain color?   2. We could hypothetically assume that the ball size follows a known distribution (in this case, zipfian). Would this help?   3. Instead of a number, how can I obtain a probability distribution of the estimated population size and volume? * * * Some clarifications over the original problem:   1. The population size is known _à priori_.   2. The population volume is also known _à priori_.   3. But the number of colors is unknown. Or the total volume of each color. Or both ratios.   4. Each sample take balls without replacement.   5. Samples are independent (after a sample is taken, all balls are replaced).


title: Chi-squared test with 0 expected values

text: My contingency table:                        heterozygous homozygous.minor homozygous.major     observed    2                 0               3     expected    0                 0               5      The expected population is composed of only AA genotype, but in the observed population we observe 2 AB genotypes. To calculate the Chi-sq for this would I just ignore the two cases where the expected = 0? So I would do: $(3-5)^2/5=0.8$


title: Ties in the data

text: Are there any "good" ways to eliminate the effects of ties in a data set for statistical analysis? For example, in finding the correlation between two series which has lots of ties? I would like to estimate the dependence between two series by a Pearson correlation coefficient and other correlation measures, like Kendall's tau. However, I don't know whether Pearson's correlation coefficient, or any other coefficient, is good for such data. The data set is around 40K observations, both variables take values in [0,1], and there are lots of ties (around 20K) at 0 and 1.


title: Evaluate text rarity in document set

text: I would like to evaluate **the rarity of each sentence in a document set**. Please let me know the state-of-the-art or a survey paper on this task. I've already checked several papers in **the document summarization problem** (e.g. [1]) which seems to be closely related my task and in fact one of the summarization methods offers me modest estimation result. But, I believe the evaluation of the text rarity should be one of the important topics in NLP field and there are some good materials to learn the topic. [1]Automatic summarising: the state of the art, Karen Sparck Jones, 2007


title: How would you explain the difference between correlation and covariance?

text: Following up on this question, How would you explain covariance to someone who understands only the mean?, which addresses the issue of explaining covariance to a lay person, brought up a similar question in my mind. How would one explain to a statistics neophyte the difference between _covariance_ and _correlation_? It seems that both refer to the change in one variable linked back to another variable. Similar to the referred-to question, a lack of formulae would be preferable.


title: Maximum Likelihood Estimator of rate parameter of the exponential distribution (MLE)

text: I want to find the maximum likelihood estimator of the "rate parameter theta of the Exponential Distribution". So i followed the following commands in `R`:                x=rexp(500,rate=2)      f <- function(x,theta){             sum(-dexp(x,rate=theta,log=T))      }      optimize(f=f,x=x,interval=c(0,5))      ### I do not understand the logic - _why do we consider`maximum=FALSE(default)` in the `optimize` function?_ ### Why haven't we used `maximum=TRUE` to find the MLE?


title: Neural network packages which allow shared weights and parallel training

text: I'm curious if there are any neural network packages out there that easily allow one to construct feed forward neural networks with shared weights, but also allow for the training to be done in parallel. Torch7 allows for easy construction of shared weights, although the parallel training support is either not there (or not documented well enough to make it obvious that it is). If it interfaced to Python, that would be even better, but this is not a requirement.


title: Is there a (parametric) distribution over sparse symmetric, binary matrix-valued random variables?

text: The title says it all, really. I'm aware of the Wishart distribution for symmetric, nonnegative-definite matrix-valued random variables, and am looking for something along these lines, but for sparse symmetric, binary matrix- valued random variables. If necessary, I'm willing to forego the requirement that the matrices be sparse. The context for this is the study of adjacency matrices for undirected graphs. In particular, “distant” vertices are unlikely to be connected, hence the sparsity (whereby “distant”, I mean dissimilarity between whatever it is that the vertices represent, measured in some suitable way). Thanks in advance.


title: Generate three correlated uniformly random variables

text: Suppose we have $$X_1 \sim \textrm{unif}(n,0,1),$$ $$X_2 \sim \textrm{unif}(n,0,1),$$ where $\textrm{unif}(n,0,1)$ is uniform random sample of size n, and $$Y=X_1,$$ $$Z = 0.4 X_1 + \sqrt{1 - 0.4}X_2.$$ Then the correlation between $Y$ and $Z$ is $0.4$. How can I extend this to three variables: $X_1$, $X_2$, $X_3$?


title: Intraclass correlation (ICC) for an interaction?

text: Suppose I have some measurement for each subject at each site. Two variables, subject and site, are of interest in terms of computing intraclass correlation (ICC) values. Typically I would use function `lmer` from R package `lme4`, and run               lmer(measurement ~ 1 + (1 | subject) + (1 | site), mydata)      The ICC values can be obtained from the variances for the random effects in the above model. However, I recently read a paper that really puzzles me. Using the above example, the authors calculated three ICC values in the paper with function lme from nlme package: one for subject, one for site, and one for the interaction of subject and site. No further details were given in the paper. I'm confused from the following two perspectives:   1. How to calculate the ICC values with lme? I don't know how to specify those three random effects (subject, site, and their interaction) in lme.   2. Is it really meaningful to consider the ICC for the interaction of subject and site? From modeling or theoretical perspective, you can calculate it, but conceptually I have trouble interpreting such an interaction.


title: What is the difference between "testing of hypothesis" and "test of significance"?

text: Is there a difference between the phrases "testing of hypothesis" and "test of significance" or are they the same? After a detailed answer from @Micheal Lew, I have one confusion that nowadays hypothesis (e.g., t-test to test mean) are example of either "significance testing" or "hypothesis testing"? Or it is a combination of both? How you would differentiate them with simple example?


title: Estimating variance for identically non independent data

text: Let $X_{ij}$ with $1\leq i<j\leq n$ (that are $X_{12},\dots, X_{1n},\dots,X_{(n-1)n}$) be ${n \choose 2}$ identically normal distributed $N(0,\sigma^2)$ such that $ \text{corr}(X_{ij},X_{rs})=\rho $ if $|\\{i,j\\}\cap\\{r,s\\}|=1$ and $0$ if $|\\{i,j\\}\cap\\{r,s\\}|=0$. I'd like to estimate $\sigma^2$ if in the large sample setting. I have proved that $\hat\sigma^2=\frac1{{n\choose 2}}\sum_{i=1}^n\sum_{j=i+1}^nX_{ij}^2$ converges to $\sigma^2$ in probability. But, how can I know that this estimator is the best estimator that we can do? For example, in the efficiency setting, do we have $\sqrt{{n\choose 2}}(\hat\sigma^2-\sigma^2)\to N(0,\tau^2)$ for some $\tau$? Any help will be appreciated.


title: Finding an equation with many variables to fit a set of data

text: I am writing a program which takes notes from a keyboard as the input, (just numbers, 1 to 88) and decides which notes are played by which hand. There are a lot of variables, for example, the position of the note(the higher it is the more chance it will be played by a right hand and vice versa), the amount of notes played by the same hand next to the target note, etc. If I have some experimental data and analysed it, what is the way to fit an equation ( a curve) to it? The most I've done of the sort is fitting a line to a set of data, so I am a bit lost. Also, as I understand, if there are two (and there probably will be more) variables than the plot is 3-dimentional , so is there still a way to do it if I have a lot of variables. I am a bit of a noob in statistics so links with detailed description etc will be very much appreciated. Thanks. EDIT: Ok, here is a little more information. For now let's take two variables, the height of the note and the number of notes played by the same hand near the given note (a radius of one bar). If the note is number 40 (Middle C) or higher it is more likely to be played by the right hand, and if the note is below 40 it is more likely to be played by the left one. I think a sigmoid function with a centre at 40 would describe it pretty well, with the shape being slightly different depending on the second variable. I have a score of around 1000 notes, each note has information about the hand it is played by, its height and the number of notes played by the same hand in the radius of 1 bar. ANOTHER EDIT: OK, here is a little more information about the nature of data I have and the purpose of all this: The purpose of the program is not analytical, but practical. A player plays a MIDI keyboard and records a MIDI file. This is basically a file containing digital descriptions of every note he played - its height, time at which it was played, length, and velocity with which it was hit. I have no information which note was played by which hand. However, I can have a good guess using the information I do know about the notes. All I need is formula describing probability of each note being played by either a left hand or a right one. The note is played by either left of right hand, so the output I am after is Left or Right, it is a bit like a coin flip problem, only this one depends on many variables rather than just 50/50. Please let me know if you need any more information. Thanks.


title: Estimating percentages as the dependent variable in regression

text: I have the rank percentages of students in 38 exams as the dependent variable in my study. A rank percentage is calculated by (a student's rank / number of students in an exam). This dependent variable has nearly uniform distribution and I want to estimate the effects of some variables on the dependent variable. Which regression approach do I use?


title: Confusion related to PCA

text: I have this confusion related to PCA. Lets say I have a 100 dimensional data, then the projection into the first principal component is given by the average over all the dimensions. How is it so?


title: Evaluating the variation in the power of a known periodicity

text: Consider the following example:               t = 1/24:1/24:365;     x = cos((2*pi)/12*t)+randn(size(t));     % if you have the signal processing toolbox     [Pxx,F] = periodogram(x,rectwin(length(x)),length(x),1);     plot(F,10*log10(Pxx)); xlabel('Cycles/hour');     ylabel('dB/(Cycles/hour');      This demonstrates the dominant periodicity in the data set. However, if I have a time series for one year worth of data, similar to that shown above, is it possible to look at how a specific frequency changes through time. For example, the time series above shows the hourly variation of a given variable throughout the year. Therefore if we knew that the dominant period was driven by a diurnal cycle (i.e. once per day) then we would know that the dominant periodicity would be 1/24. So, if we know the dominant periodicity to be 1/24, is it possible to see how the power of this specific periodicity changes in time (i.e. throughout the year)?


title: Principal Components, Canonical Correlation and Eigenvalue problems

text: It is well known that the solution to the optimization problems proposed in Principal Components and Canonical Correlation Analysis are given by the solution to eigenvalue problems and generalized eigenvalue problems respectively. My linear algebra is quite rusty and I wanted to know if anyone could give me a semi-formal, albeit intuitive explanation on why this is so.


title: Two-sided permutation test vs. two one-sided

text: I am working with a set of genes for which I have both methylation values $\beta$ (continuous on unit interval) and gene expression $E$ (non-negative continuous) that I want to test for correlation $\hat\rho(log_2E,\beta)=R$. To filter out only the significant correlations I perform a permutation test where for each site I randomly permute $E$ a large number of times and calculate $R$. This gives me a distribution under the null hypothesis $p(R|H_0)$, i.e. when there is no dependence between $\beta$ and $E$, that I can compare $R$ against. Now to my question, I want to test both for positive and negative correlation. Do I make two tests per gene and use double multiple testing correction, where $p_- = \int\limits_{-\inf}^Rp(x|H_0)dx$ $p_+ = \int\limits_{R}^\inf(x|H_0)dx$ or can and should I calculate a two-tailed p-value directly based on $|R|$ or such, and use normal multiple testing correction? I can only think of the following, but it feels wrong and has lousy power. $p_2 = \int\limits_{-\inf}^{-|R|}p(x|H_0)dx + \int\limits_{|R|}^\inf(x|H_0)dx$ ![enter image description here](http://i.stack.imgur.com/nD8CJ.png) **Edit:** Updated the figure and added the missing x-axis label. The y-axis is the same for all panels.


title: ADF Test on Volatile Series. Expected results?

text: I tried performing the Augmented Dickey-Fuller Test on SP500 data from MASS library and got the following results:               adf.test(SP500)          Augmented Dickey-Fuller Test          data: SP500     Dickey-Fuller = -13.5835, Lag order = 14, p-value = 0.01     alternative hypothesis: stationary      The results seem to indicate that the SP500 series is stationary. However by definition stationary imply constant mean and variance. My question is how can the test conclude that the series is stationary when the series exhibit blatant volatility which imply that the variance of the series is not constant? Also can someone tell me what the parameter alternative and k mean here:               adf.test(x, alternative = c("stationary", "explosive"), k = trunc((length(x)-1)^(1/3)))      I am still stuck on this question, are financial returns time series supposed to be stationary?because the ADF test is telling so but I can't seem to agree.