title: Using MCMC to evaluate the expected value of a high-dimensional function

text: I am working on a research project that is related to optimization and recently had an idea to use MCMC in this setting. Unfortunately, I am fairly new to MCMC methods so I had several questions. I'll start by describing the problem and then asking my questions. Our problem boils down to estimating the expected value of a cost function $c(\omega)$ where $\omega = (\omega_1,\omega_2,...\omega_h)$ is an $h$-dimentional random variable with a density $f(\omega)$. In our case, a closed form version of $c(\omega)$ does not exist. This means that we have to use Monte Carlo methods to approximate the expected value. Unfortunately, it turns out that estimates of $E[c(\omega)]$ that are generated using MC or QMC methods have too much variance to be useful within in a practical setting. One idea that we had to use an importance sampling distribution to generate sample points that will produce a low variance estimate of $E[c(\omega)]$. In our case, the ideal importance sampling distribution, $g(\omega)$, has to be roughly proportional to $c(\omega)f(\omega)$. Seeing how $g(\omega)$ is known up to constant, I am wondering whether I can use MCMC along with the proposal distribution $c(\omega)f(\omega)$ to eventually generate samples from $g(\omega)$ . My questions here are:   * Can MCMC be used within this setting? If so, what MCMC method would be appropriate? I am working in MATLAB, so I have a preference to anything that already has a MATLAB implementation.    * Are there any techniques that I can use to speed up the burn-in period for MCMC. And how can I tell that the stationary distribution has been reached? In this case, it actually takes a fair bit of time to calculate $c(\omega)$ for a given $\omega$.


title: Non-seasonal periodicity in Time Series residuals

text: I have been working on a forecast model in Excel extrapolating from a small (150 data points) monthly time series. I've converted into a year/year percentage change series to get it stationary, deseasonalized it with monthly dummies, and accounted for serial-correlation (it has a two period error correlation). However, after graphing the residuals like so... ![Graph of actual vs. predicted errors](http://i.stack.imgur.com/HbDdH.png) ...it's very clear that there is non-seasonal, semi-regular periodicity that still needs to be removed. What is the recommended approach for modeling this kind of error (hopefully one that doesn't require R)? Edit: The PACF for the first two lags were (1) -0.481, and (2) -0.2786. PACF values for lags greater than were not significant.


title: What is predictor offset in linear regression?

text: I've been reading the post about removing intercept that boosts $R^2$ (HERE). An answer to that post said that: > In essence, that means our predictor had better have a strong mean offset > itself and that this mean offset should dominate the variation of the > predictor. My question is: what is "a predictor mean offset"? I've searched the Internet and found nothing :(


title: Proper "normalization" of post-test data to pre-test conditions

text: I'm performing an _in vitro_ assay in which cells are treated with different combinations of two different compounds and the activity of the cells following treatment is measured. There are 3 concentrations of compound A and 4 levels of compound B (one concentration for each compound is 0, as a control), and all permutations are used. Cells have been plated in microplates, and in theory there should be an equal number of cells in every well. Four wells receive each of one of the 12 different compound combinations. To account for preparation inconsistencies (i.e., unequal numbers of cells between wells) before treatment, the same activity measurement method is done before treatment, since activity is proportional to cell number. The question I want to answer is: how do these compounds increase (or decrease) cell activity? The pre-test measurement did find differences in the average activity between groups (wells were assigned to groups via the plate layout before adding the cells). I want to account for this difference prior to treatment so that treatment effects are not clouded by well-to-well variations. I've read that an ANCOVA is a good way to deal with this. I wonder if other normalization methods would work. I considered determining the maximum activity level of all the groups before treatment, dividing each group's pre- treatment average by that maximum to create a normalization factor, and then dividing the post-test measurements by the normalization factor. Is this inappropriate, and what would be a good way of accomplishing my task?


title: Counts of binned data by group

text: I have continuous data $A$ and categorical data $O$. I need counts of $A$ in bins by group $O$. I'm working in R. I know how to bin data (using `cut2`) and how to get the counts of $O$ (using `aggregate` or `by` or `describe` or `summaryBy`). I could also get what I need by running one of these functions with a subset (bin) of the data one by one, but is there a way to do this automatically in one command? **Update:**   The output I want is just basic descriptive statistics (counts in this case, but also means and/or medians for example):               Age (A)................Outcome (O)     ........................Yes....No     18-27 ...................10....12     28-37 ....................2....11     38-47 ....................9.....7      For example:               by(subset(X, X$Age>47)$Age, subset(X, X$Age>47)$Outcome, length)      gives me the number of data points for ages > 47 where `Outcome` is true and where `O` is false. I want to do this for all bins automatically, without having to run separate commands with long > and < conditions for each bin.


title: odds ratio: the purpose and interpretation

text: Could someone please explain the purpose of the odds ratio and how it could be interpreted (i.e., the origin of their usage is in primary interest for me)? Why don't people simply use the difference between two proportions instead of odds ratio? Reference to the proper mathematical book will be great.


title: Import stock price from Yahoo Finance into R?

text: I would like to import the "Last Trade" stock price from Yahoo finance into R. The intention is to work with (almost) real time data. Are there any solutions? Thanks in advance for any helpful comment.


title: Why is it that a lower R-Squared on a difference regression model could be better than higher R-squared on a levels regression model

text: I am reviewing a time series regression model that uses the log of the year over year change in sales as the dependent variable and the log of the year over year change in another economic index as the independent variable. The regression model produced an R-squared of .60. When I expressed concern about the low R-squared, I was told it was a very good R-squared given the use of log differences instead of levels and that I can't really compare it to an R-squared derived from using levels. I thought R-squared represented goodness of fit. So why would a regression with a 0.60 R-Squared based on difference be better than a regression with a 0.85 R-sqared based on levels. Thanks for your input.


title: Should one remove highly correlated variables before doing PCA?

text: I'm reading a paper where author discards several variables due to high correlation to other variables before doing PCA. The total number of variables is around 20. Does this give any benefits? It looks like an overhead to me as PCA should handle this automatically.


title: Smoothing objective function

text: I have a current problem where i have a continuous error real function of n-variables, which is computationally expensive and (possibly) has multiple local minima (optima). I want to somehow smooth the error surface and decrease the number of local minima with as low deformation as possible (may be controlled by a parameter) to the surface so that the global minima and main features of the surface is preserved. After which I can apply gradient decent or simulated annealing or other algorithms on it to get a better algorithm. Can anyone point where to start on smoothing so that i can try to achieve what i want? Note that I am not knowledgeable enough in optimization techniques and statistics, but do point me to any resources/poitners which i can study to get into it.


title: Multiplicative mixed models for analyzing variety-by-environment data

text: I got some nice data on multi-environmet trials (MET) data for genotype evaluation and would like to use some new developed techniques as discussed in Smith et. al 2005 Smith et. al 2005. I'm specifically interested in **Factor Analytic (FA) structure**. Authors mentioned codes for these methods will be available on request. > There are now several statistical packages (including ASReml, GENSTAT, > S-language packages and SAS; Littel et al. 1996) that allow REML estimation > of a range of mixed models. The present authors have found the packages > ASReml and GENSTAT and the samm functions (through S-language environments) > to be the most suitable for the analysis of MET data, both in terms of the > generality of models that can be fitted and the ease with which predictions > and inference about varietal effects can be made. All models in the current > paper are easily fitted and summarized using these software (code is > available from the authors on request). Even after multiple requests to the authors I haven't heard back from the authors. I'm wondering if someone has tried these models and be kind enough to share worked examples, that would be great. I'm looking forward for some positive response. Smith et al. (2001b) used the following mixed model version of multiplicative models ${\small \begin{eqnarray*} (\boldsymbol{I}_{e}\otimes\boldsymbol{I}_{g})\boldsymbol{\eta}_{eg\times1} & = & (\boldsymbol{1}_{e}\otimes\boldsymbol{1}_{g})\mu+(\boldsymbol{I}_{e}\otimes\boldsymbol{1}_{g})\boldsymbol{E}_{e\times1}+(\boldsymbol{1}_{e}\otimes\boldsymbol{I}_{g})\boldsymbol{G}_{g\times1}+\underbrace{(\boldsymbol{I}_{e}\otimes\boldsymbol{I}_{g})(\boldsymbol{GE})_{eg\times1}}\\\ & = & (\boldsymbol{1}_{e}\otimes\boldsymbol{1}_{g})\mu+(\boldsymbol{I}_{e}\otimes\boldsymbol{1}_{g})\boldsymbol{E}_{e\times1}+(\boldsymbol{1}_{e}\otimes\boldsymbol{I}_{g})\boldsymbol{G}_{g\times1}\\\ & & +\underbrace{(\boldsymbol{\Lambda}_{E_{e\times k}}\otimes\boldsymbol{I}_{g})\:\boldsymbol{f}_{G_{kg\times1}}+(\boldsymbol{I}_{e}\otimes\boldsymbol{I}_{g})\delta_{eg\times1}}\end{eqnarray*} }{\tiny }$ where $\boldsymbol{\Lambda}_{E_{e\times k}}$ is a matrix of environment loadings, $f_{G_{kg\times1}}$ is the associated vector of genotype scores and _k_ is the number of components (multiplicative terms) included in the model. The authors assumed that the environments are fixed and $\boldsymbol{G}_{g\times1}$, $\boldsymbol{f}_{G_{kg\times1}}$, and $\delta_{eg\times1}$ are random effects with $\small \left(\begin{array}{l} \boldsymbol{G}_{g\times1}\\\ \boldsymbol{f}_{G_{kg\times1}}\\\ \boldsymbol{\delta}_{eg\times1}\end{array}\right)\sim\mathcal{N}\left(\begin{array}{ccc} \left[\begin{array}{c} \boldsymbol{0}\\\ \boldsymbol{0}\\\ \boldsymbol{0}\end{array}\right] & , & \left[\begin{array}{ccc} \sigma_{g}^{2}\,\boldsymbol{I}_{g} & \boldsymbol{0} & \boldsymbol{0}\\\ \boldsymbol{0} & \boldsymbol{I}_{k}\otimes\boldsymbol{I}_{g} & \boldsymbol{0}\\\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{\Psi}_{e}\otimes\boldsymbol{I}_{g}\end{array}\right]\end{array}\right)$ $\small E((\boldsymbol{I}_{e}\otimes\boldsymbol{I}_{g})\boldsymbol{\eta}_{eg\times1})=(\boldsymbol{1}_{e}\otimes\boldsymbol{1}_{g})\mu+(\boldsymbol{I}_{e}\otimes\boldsymbol{1}_{g})\boldsymbol{E}_{e\times1}$ ${\small \begin{eqnarray*} \mathrm{var}((\boldsymbol{I}_{e}\otimes\boldsymbol{I}_{g})\boldsymbol{\eta}_{eg\times1}) & = & \sigma_{g}^{2}\,(\boldsymbol{J}_{e}\otimes\boldsymbol{I}_{g})+(\boldsymbol{\Lambda}_{E_{e\times k}}\otimes\boldsymbol{I}_{g})(\boldsymbol{I}_{k}\otimes\boldsymbol{I}_{g})(\boldsymbol{\Lambda}_{E_{e\times k}}\otimes\boldsymbol{I}_{g})^{T}+\boldsymbol{\Psi}_{e}\otimes\boldsymbol{I}_{g}\\\ & = & \sigma_{g}^{2}\,(\boldsymbol{J}_{e}\otimes\boldsymbol{I}_{g})+(\boldsymbol{\Lambda}_{E_{e\times k}}\boldsymbol{\Lambda}_{E_{e\times k}}^{T}\otimes\boldsymbol{I}_{g})+\boldsymbol{\Psi}_{e}\otimes\boldsymbol{I}_{g}\\\ & = & \sigma_{g}^{2}\,(\boldsymbol{J}_{e}\otimes\boldsymbol{I}_{g})+(\boldsymbol{\Lambda}_{E_{e\times k}}\boldsymbol{\Lambda}_{E_{e\times k}}^{T}+\boldsymbol{\Psi}_{e})\otimes\boldsymbol{I}_{g}\\\ & = & (\underbrace{\sigma_{g}^{2}\,\boldsymbol{J}_{e}+\boldsymbol{\Lambda}_{E_{e\times k}}\boldsymbol{\Lambda}_{E_{e\times k}}^{T}}+\boldsymbol{\Psi}_{e})\otimes\boldsymbol{I}_{g}\\\ & = & (\underbrace{\boldsymbol{\Lambda}_{E}^{*}\boldsymbol{\Lambda}_{E}^{*T}}+\boldsymbol{\Psi}_{e})\otimes\boldsymbol{I}_{g}\end{eqnarray*} }$ where $\boldsymbol{\Psi}_{e}$ is a diagonal $e\times e$ matrix with elements commonly referred to as specific variances. Thanks


title: Analysis of rates and post hoc test with offset variable using glmer.nb

text: I am investigating variation in pollinator visitation rate (number of visits per inflorescence) with treatment and time category as fixed factors. Block is a random factor. Following Zuur et al. (2009), I used the number of visits as the response variable with the $\log$(number of inflorescences) as the offset variable in the analysis. A Poisson model was overdispersed, and therefore I opted for a negative binomial model in `lme4`, as follows:   `model1 = glmer.nb(visits ~ treat + timecat + offset(log(infl)) + (1|block))` I am specifically interested in differences in visitation rates between treatments. I therefore performed a post hoc test: `OPexp1 = glht(model1,mcp(treat = "Tukey")); plot(cld(OPexp1))` When I plot these results, I get number of visits on the $y$ axis. But what I want is visitation rate (visits per inflorescence). Is the post hoc test taking into account the number of inflorescences? Or how do I specify that the post hoc test should be performed using visitation rate? I thought that running the model above is analogous to analysing the visitation rate, so are the Tukey's differences between treatments for number of visits also the same for visitation rates? I assume what is happening is that fitted values are currently expressed as $μ × V$ (where $μ =$ visitation rate and $V =$ number of inflorescences), but how do I specify that they should be expressed as $μ$ (visits per inflorescence) only? On p. 240, Zuur et al. (2009) mentions that this is possible, but I have not been able to find an example. I am quite new to R, and this is my first post here, as I have been unsuccessful in finding an answer elsewhere, so any advice or kick in the right direction would be much appreciated. Kind regards.


title: Confidence Interval Coverage-error and Type I error

text: Could somebody explain to me the relationship between coverage error and type one errors in multiple comparisons testing, if there is one in fact? Does a coverage error occur when the true value of the comparison is in fact not within the parameters of the confidence interval? Further, does this have any influence on the probability of a type I/II error for that comparison? Thanks


title: What is the number of free parameters for a directed acyclic graph?

text: I am trying to use Maximum Likelihood Estimation to learn the structure of a DAG, G. How is the number of free parameters of G calculated to compare the complexity of different graphical models? Is it based on one of the following, or something else?   * Number of edges in the graph   * Maximum number of possible edges in the graph   * Maximum number of parents (children)


title: Probability that the probability for a binomial distribution is below a certain value

text: This is a slightly odd problem I ran across recently. Assuming a population of size $N$ following a binomial distribution with unknown $p$, how many different samples must be taken from the population to say that $P[p<\phi]>\alpha$, assuming all the samples taken turn out to be failures? My approach was to state that the probability of $p=x$ given $n$ failures is $$ P[p=x|n] \propto (1-x)^n $$ Normalizing this so that $\int_0^1 P[p=x|n]\ dx=1$ gives $$ P[p=x|n] = (n+1)(1-x)^n $$ The probability that $p<\phi$ is then easy to calculate as $$ P[p<\phi|n] = \int_0^\phi (n+1)(1-x)^n\ dx = 1 - (1 - \phi)^{n + 1} $$ And finally solve the inequality $P[p<\phi|n]>\alpha$ $$ 1 - (1 - \phi)^{n+1} > \alpha \ \Rightarrow\ n > \frac{\ln 1 - \alpha}{\ln 1 - \phi} - 1 $$ Now, my first question regarding this is whether my approach is correct or not. Secondly, is there an intuitive way to explain/understand the equations above? Above I solved it only mathematically, but I'm really more interested in the _why_ than the how. For example, $1-(1-\phi)^{n+1}$ above, given $n$ failures, reads to me like the probability that the next sample taken is _not_ a failure. But I see no connection between that and the left hand side $P[p<\phi|n]$.


title: Averages of random subsets of variables

text: Consider data where each observation was generated as follows.   * We draw $Z_1,...,Z_m$ from some distribution. (Possibly they're independent or related in some other simple way.)   * Next, based on the $Z_1,...,Z_m$, we choose a sequence $0=I_0 < I_1 < ... < I_N=m$ so that, for each $k$, (i) $I_k-I_{k-1}$ is not too small and (ii) the sample variance within $Z_{I_{k-1}+1},...,Z_{I_k}$ is small. (I am intentionally somewhat vague here - I am open to making various different assumptions along these lines.)   * We generate the observed variables $X_1,...,X_N$ as $X_k=$ the average of $Z_{I_{k-1}+1},...,Z_{I_k}$. For example, the hidden sequence $Z=(0.1, 0.3, 0.2, 1.3, 1.2, 0.1)$ might lead to the observed sequence $X=(0.2, 1.25, 0.1)$ [or perhaps to $X=(0.2,0.86)$ due to (i) above]. Does anyone here know whether this type of setup has been studied before, and if so, what are some keywords to search for or papers/books to look at? Thanks in advance for any answers! Added on Apr 21: The motivation is as follows. Think of each $Z$ sequence as SNP data from a single patient. In order to anonymize the data for a public release a procedure like the one I described above can be performed. Based on the anonymized data $X$, I want to predict survival and/or to identify SNPs that are relevant to survival. Note that $I$, $N$, and $X$ are all functions of $Z$, so they will be different for each patient. Also note that that $I$'s are observed, i.e., I know which $Z_j$'s were averaged to produce each $X_k$.


title: Frequency on Tbats function in R

text: I have 3 complete years + 4 weeks of weekly time series data, of which one of the years is a leap year. To calculate its frequency should I do (365x2+366)/(3*7)? Thanks!


title: Why are these MLE estimates biased?

text: I estimate the parameters of survival data with censoring which is simulated from Weibull distribution. The mean time to event was set to 10 by choosing the combinations of shape and scale parameters as follows: (constant hazard,shape=1 scale=10) and (decreasing hazard,shape=0.6 scale=6,646). The censoring times were simulated from exponential distribution.In order to attain target censoring rate(0.76),different mean of exponential distribution was set. For constant hazard, the mean of exponential was set to 3.1, and 1 for decreasing hazard. The mean time was calculated when the parameters were obtained in R by `survreg()`: mean of weibull= scale* gamma(1+1/shape). After 1000 simulations, it was supposed that the mean of mean time for constant hazard is not far from the one for decreasing hazard, because their censoring rate is the same(0.76). However, it surprises me that the mean of mean time for constant hazard is 10.9, and the one for decreasing hazard is 15.1. So I am wondering that if the constant hazard make the mle estimates more biased than decreasing hazard? The R code is below:               library(survival)     n<-100            # the number of patients      foltime<-12       # calendar time when study ends     simulations<-1000     mean_time_simulation<-numeric()     censored_proportion<-numeric()     mean_time_complete<-numeric()     shape_input<- 0.6   #or shape=1,scale=10     scale_input<-6.646  #or shape=1,scale=10     mean_time <- scale_input * gamma(1+1/shape_input)       mean_time_censored <- 1       for(counter in 1:simulations){     te<-rweibull(n,shape_input,scale_input)      # patient times Weibull     tc<-rexp(n, rate = 1/mean_time_censored )   # times when patients censored     t0<-runif(n,0,3)      #study times start     t_end<-t0+te          # calendar time patient death     t_calendar_censored<-t0+tc      #  calendar time when patient censored     it<-c(1*(t_end<t_calendar_censored)*(t_end<foltime))    #  0 = censored, 1 #dead     time1<-(pmin(t_calendar_censored, rep(foltime,n))-t0)*    (t_end>pmin(t_calendar_censored,rep(foltime,n)))    +   te*(t_end<=pmin(t_calendar_censored,rep(foltime,n)))     censored_proportion[counter] <- (n - sum(it)) / n  # proportion of patients censored     model_complete<-survreg(Surv(time1,it)~1,dist="weibull")     lambda_complete<-exp(summary(model_complete)$table[1])         gamma_complete<-1/exp(summary(model_complete)$table[2])     mean_time_complete[counter]<-lambda_complete*gamma(1+1/gamma_complete) #mean survival  time     bias_mean_time_complete<- mean(mean_time_complete) - mean_time     }     mean(mean_time_complete)


What is validity of an inference procedure? The inference procedure can be point estimator, confidence interval, or a testing rule. I saw "validity" in Diggle et al's Longitudinal Data Analysis: ![enter image description here](http://i.stack.imgur.com/PF1uH.png) Thanks!


title: Understanding multinomial distribution

text: There are K categories. Let x be a discrete random variable taking on values $1,2...,K$. Given data $\textbf{X}= \\{x_1, x_2, ..., x_N\\}$, and suppose $\textbf{p}$ is a vector probabilities where, $p(x= k) = p_k$. Normally, we see the following:   1. Suppose $N_k = \sum_i I(x_i = k)$, i.e. count number of times outcome k occurs in n trials. And we model the set of counts $N_k$ as multinomial distribution.  $$p(N_1...N_k | \textbf{p}) = \frac{N!}{N_1!...N_K!} \prod^K_{k=1} {p_k}^{N_k}$$ This is what we see normally. However, I came across something interesting on the web and my question is the following: can the following be a multinomial distribution ? And why ? What if we model $\textbf{X}$ as multinomial distribution. $$p(\textbf{X} | \textbf{p}) = \prod^{K}_{k=1} {p_k}^{N_k}$$ Is it still a valid pmf for multinomial distribution ? At first ,I thought it was impossible because it lacks permutation, which is the normalizing constant that makes the pmf sum to 1. But the research paper indicates both of these are multinomial distribution, although they are modelling different things.


title: When using the exreg argument in auto.arima in R, is the main series still being used as a predictor?

text: What is the purpose of and how to use the xreg argument when fitting ARIMA model in R? My question relates to how the xreg argument is working when fitting an ARIMA model (purpose of xreg described in the post above.)               library(forecast)     fity <- auto.arima(y, xreg=x)      In the code above, is the ARIMA model being fit using past values of $y$ as predictors as well as the past values of $x$? i.e. is it assumed that the past $y$ values would be used as they would be if running an ARIMA modelling without an explanatory variable $x$? Or if I use $\text{xreg} = x$ and want the $y$ series to be used as predictors, would I have to add $y$ to the xreg matrix to force that to occur? Thanks.


title: Correlated/dependent random variables

text: This may be a very basic question or a very difficult one, I am not sure. What does it mean if two random variables are perfectly correlated? What does it mean if two random variables are perfectly dependent (whatever that means)? When two random variables are independent joint pdf can be factorized into marginal pdfs. But if they are say perfectly correlated, what does that say? Does it mean they come from distributions with same functional forms but with different variances? What does it mean to be perfectly dependent? Conditional distribution collapses to a dirac delta function? Thanks


title: Connection between sum of normally distributed random variables and mixture of normal distributions

text: If you have two independent random variables that are normally distributed (not necessarily jointly so), then their _sum is also normally distributed_ , which e.g. means that its excess kurtosis is $0$. On the other hand in case the mixture of one-dimensional normal distributions the mixture distribution can display non-trivial higher-order moments such as skewness and kurtosis (fat tails) and multi-modality, even in the absence of such features within the components themselves (see also this video for an easy enough example). This means that the result _doesn't have to be normally distributed_. **My question**   How can you reconcile these two results and what is their connection?


Is Baron and Kenny method for mediation now outdated? This is because I got a comment from reviewer saying asking me to refer to Rucker, D.D., Preacher, K.J., Tormala, Z.L. & Petty, R.E. (2011). Mediation analysis in social psychology: Current practices and new recommendations.


title: Massive Textual Corpus

text: > **Possible Duplicate:**   >  Where to find a large text corpus? I know someone has asked a similar question here, but I'm wondering whether anyone knows of a large textual corpus that is available for research use. The number of documents it contains isn't terribly important--rather, I'm looking for something that is on the TB-/PB-level that I can use to test the scalability of some algorithms. I thought about using the english Wikipedia data dump, but I think it's only about 25GB. My other thought was using a database of twitter messages, but, if I remember right, those aren't freely- available. Does anyone have a recommendation?


title: How to get from joint distribution F(x,y) to f(x,y) to calculate the marginal distribution of X?

text: I have this homework question I'm not 100% sure how to tackle. I have a random vector with joint distribution function F(x,y) and am asked to find the marginal distribution function. I think need to get to $fx(x,y)$, but I'm not sure how to do this. My thinking is that I need to differentiate F(x,y) by both x and y, to get f(x,y) and then integrate over y, but this gets me the marginal density for x. Any pointers would be greatly appreciated EDIT - Thanks Glen. I've switched to my proper user id. It's late here, I'll revisit tomorrow. EDIT2 - to clarify, the first q in the homework starts at $F(x,y)$ and asks for marginal _probability_. The second starts with $f(x,y)$ and asks for marginal _density_ so I'm assuming that in the first instance, they're after a pdf, and in the second, a cdf; to rule out suspicions of casual terminology use.


title: Is a sample covariance matrix always symmetric and positive definite?

text: When computing the covariance matrix of a sample, is one then guaranteed to get a symmetric and positive-definite matrix? Currently my problem has a sample of 4600 observation vectors and 24 dimensions.


title: Concerning WN, IID tests

text: If we have a null hypothesis and subject it to Turning Point test, rank, runs portmant. etc. If the p value fall below 0.05 in some cases and not for others, do we then say that we reject the hypothesis? I.e. what can we say about the WN hypothesis if our series passes some tests but fails at others?


title: What is the MLE if the variance is heteroskedastic

text: The normal distribution of MLE is start ${x_1,...,x_n}$ random sample from $\sim\mathcal{N}(\mu,\sigma^2)$ This family of distributions has two parameters: $\theta = (\mu, \sigma)$, so we maximize the likelihood over both parameters simultaneously using the logarithm of the likelihood. I know this. But what if ${x_1,...,x_n}$ is a random sample from the normal distribution $\sim\mathcal{N}(\mu, \sigma_i^2)$ ? This means the variance is not homoscedastic. Than what is the MLE of $(\mu, \sigma_i^2) $?


title: Adding decays to time series in R

text: Working with my time series data I need to add decays to a certain variables. That is, when the variable consists of a series of impulses, I need to prolong these effects of these in order to better depict the true underlying nature of the effect. I would also like to experiment with different types of decay lengths (number of observations before cut-off) and types (e.g. linear and exponential). Moreover, the code should be able to cope with summing decays when impulses are within each others decay length. Demo data               time  impulse      1      0      2    100      3      0      4      0      5      0      6      0      7      0      8    200      9    100     10      0     11      0     12      0     13      0     14      0     15      0     16      0     17      0     18     60     19      0     20      0


title: Statistical tests to do data quality checks

text: I want to identify bad responses in the data we get from our vendors. To assess the quality of the data, what statistical tests can i apply to identify bad respondents? I performed Discriminant analysis to identify the % similarity of responses among surveys completed by a particular interviewer. How can i further authenticate my findings? (added from comments) What i mean by bad responses is: Interviewers who have very high similarity in the responses they collect in the survey. they might tick for "4" for a particular question for all the respondents (Flat Liners or Speeders) Suppose an interviewer has conducted 10 surveys and i notice similarity in the answers of those questionnaires for each question. Is there any statistical test which authenticates my finding that the interviewer has probably filled responses himself or has ticked the same response for each question? How can we check a particular data using Statistical tests? I hope i have clarified my question.


title: Getting a time plot on levels using the first differences of a time series in R

text: In my current regression model, one of the explanatory variables I need is "capital" (K). The country under consideration here is Côte d'Ivoire (West Africa) and, as you may already know, finding data for small open economies is very often an extremely difficult task. I was able to find the data for "capital formation", which is essentially the first differences of the actual capital values (dK), but I am required to show the plot on levels (and not only on differences). It is obvious that it is impossible to find the actual values from just the first differences, but there may be a way to have an exact similar plot using the differences as a basis. How should I go about it using R? Your answer will be much appreciated.


title: Mean versus imputation for missing data in the case of an ordinal scale

text: Is mean or mode better for replacing missing data for an ordinal scale? I'm thinking mode is better because the respondent has to choose between integer values (1, 2 and so on) bu I am wondering is mode imputation doesn't create bias by favoring the value that appears most often. Thank you!


title: Standard Error In Logistic Regression

text: I recently asked a friend what an "acceptable" standard error is for a coefficient in conditional logistic regression, and I got the logical, "it depends" response. Knowing that every model or study is different, what is the "rule of thumb" for minimum standard error for coefficients in logistic regression? I was once told 0.05... but haven't been able to find similar statements on the Google machine. Thanks!


title: Co-occurrence of properties in a population

text: I have 150 properties that may occur in a population of 10000 people. Individual people may have none, one or a couple of these properties. The properties are not mutually exclusive and have different frequencies in the population. I would like to answer two questions here:   1. Are particular sets of these properties associated, that is do they occur together in a person more often than can be expected by chance? The method should not only analyse for pair-wise co-occurrence of two properties, but also identify sets of properties that occur together. Ideally I'd like to have a probability measure like a p-value for a given set of properties telling me the likelihood of observing this combination of properties just by chance.   2. Given a sample of e.g. 50 out of the 10000 people I'd like to know whether the co-occurrence of properties observed in the sample is significantly different from a random sample of people from the population. How would I best address this in statistical / mathematical terms? Are there any tools you could recommend for the calculation?


title: Combining confidence intervals

text: I have a collection of efficiency curves i.e. numbers between 0 and 1 as a function of a physical variable. Each efficiency point on the curve has an associated Clopper Pearson 1-sigma confidence interval. I now need to combine these to obtain a total efficiency as a function of said variable. This means I multiply all the corresponding efficiency points together so that eff_total for point i is eff_1*eff_2*...*eff_n, etc. How do I propagate the uncertainty in the form of these confidence intervals to eff_total?


title: When do Markov random fields $\neq$ exponential families?

text: In their textbook, Graphical Models, Exponential Families and Variational Inference, M. Jordan and M. Wainwright discuss the connection between Exponential families and Markov Random Fields (undirected graphical models). I am trying to understand better the relationship between them with the following questions:   * Are all MRFs members of the Exponential families?    * Can all members from the Exponential families be represented as an MRF?   * If MRFs $\neq$ Exponential families, **what are some good examples of distributions of one type not ncluded in the other** ?  From what I understand in their textbook (Chapter 3), Jordan and Wainwright present the next argument: * * *   1. Say we have a a scalar random variable X that follows some distribution $p$, and draw $n$ i.i.d. observations $X^1, \ldots X^n$, and we want to identify $p$.   2. We compute the empirical expectations of certain functions $\phi_\alpha%$  $\hat{\mu}_\alpha= \frac{1}{n}\sum^n_{i=1}\phi_\alpha(X^i), $ for all $\alpha \in \mathcal{I}$ where each $\alpha$ in some set $\mathcal{I}$ indexes a function $\phi_\alpha: \mathcal{X} \rightarrow R$   3. Then if we force the following two sets of quantities to be consistent, i.e. to match (to identify $p$):     * The expectations $E_p[(\phi_\alpha(X)]=\int_\mathcal{X}\phi_\alpha(x)p(x)\nu(dx)$ of the sufficient statistics $\phi$ of the distribution $p$      * The expectations under the empirical distribution we get an **underdetermined problem** , in the sense that there are many distributions $p$ that are consistent with the observations. So we need a principle for choosing among them (to identify $p$). If we use the principle of maximum entropy to remove this undeterminancy, we can get a single $p$: $\DeclareMathOperator*{\argmax}{arg\,max} p^* = \argmax_{p\in{\mathcal{P}}} \,H(p)$ subject to $E_p[(\phi_\alpha(X)] = \hat{\mu}_\alpha$ for all $\alpha \in \mathcal{I}$ where this $p^*$ takes the form $p_\theta(x) \propto $ exp${\sum_{\alpha \in \mathcal{I}}\theta_\alpha \phi_\alpha(x)},$ where $\theta \in R^d$ represents a parameterization of the distribution in exponential family form. In other words, if we   1. Make the expectations of the distributions be consistent with the expectations under the empirical distribution   2. Use the principle of maximum entropy to get rid of undetermination $\rightarrow$ We end up with a a distribution of the exponential family. * * * However, this looks more like an argument to introduce exponential families, and (as far as I can understand) it does not describe the relationship between MRFs and exp. families. Am I missing anything?


title: Correcting for multiple 2-way anova testing

text: I want to study the effect of two factors on a paramater that I measure in several regions. I'm running a 2-way ANOVA test for each region (about a dozen) and I am wondering how I can deal to avoid the increase of the type I error. I'm thinking about performing an FDR correction on the p-values that I get with the ANOVA tests. I get 3 p-values for each test (one for each factor and the third for interaction between factors), so approximately 36 p-values altogether. I wonder how to appropriately correct for multiple tests and particularly whether I should correct the p-values all together or separately for each factor and for the interaction. Thank you !


title: Estimating number of data points in specific range from percentile data

text: If I have a dataset broken down into percentiles, such that I know there are _n_ values, and 10% of these values fall below _x1_ , 20% are below _x2_ and so on, how can I estimate an answer to the following question: How many values fall below _x_ , for arbitrary _x_? Or more specifically: How many values are between _x_ and _y_? Note that the intervals between percentiles are not uniform. In case I am on the wrong track, and in case it helps, here is the actual problem: * * * I have a set of data for UK gross annual earnings per job, broken down into percentiles: 10 20 25 30 40 50 60 70 75 80 90 91 92 93 94 95 96 97 98 99 I want to calculate the gross income tax for any arbitrary set of tax bands, e.g., how much revenue would the Treasury receive in one year if income tax rates were 0% on the first £10k, 10% on £10k-£20k, 15% on £20k-£30k, 25% on £30k+


title: Proper Number of Clicks for Conversion Rate Testing

text: I need to track the conversion rate of a landing page. I want to figure out how many clicks I need to get before I can be sure that I have an accurate conversion rate calculation. I first figured I would calculate the confidence interval using the appropriate formula for a proportion, and find the number of clicks needed to make the interval sufficiently tight. However, according to this link you need at least 1,000 clicks before you can assume the distribution is normal. Unfortunately I can't get that many clicks for all of my landing pages. (I'm using landing pages as a parallel for what I'm really doing - in reality I'm testing hundreds of CPA advertisements for their conversion rates which prevents me from giving each thousands of clicks. This detail doesn't change the problem, though, so it can be ignored). I then had a different idea - I'd go through the click/conversion data that I do have so far for previous landing pages using different "window sizes", such as a window of 10 clicks, 50 clicks, 100 clicks, 200 clicks, etc. For each window size, I'd go through each click in sequential order and calculate the conversion rate using the next x clicks, starting from the one I'm inspecting. I'd then calculate the variance of the conversion rates I found. I'd then find a window size that gives me a desired amount of minimal variance across the majority of the landing pages, and then use that as the number of clicks to test for for future tests. Does this approach seem to make sense / has anyone heard of something similar? Probably more importantly, what is the normal way to accomplish this goal? Thanks!


title: What is the p-value for paired t-test if the two set of data are identical

text: As title, I have two set of data that are identical and performed paired t-test using R, the p-value for the difference estimated (which is essentially 0) is NaN. I wonder if there can be a theoretical p-value for this case.


title: Simulation about quantile regression

text: This is what I have done in R:               > quantile<-rq(log(response)~gender+age+I(age^2)+gender*age+gender*I(age^2),data=Modified,     + tau=0.9)     > summary(quantile)          Call: rq(formula = log(response) ~ gender + age + I(age^2) + gender *          age + gender * I(age^2), tau = 0.9, data = Modified)          tau: [1] 0.9          Coefficients:                         Value    Std. Error t value  Pr(>|t|)     (Intercept)          3.50682  0.04854   72.24571  0.00000     genderMale           0.05580  0.05556    1.00419  0.31534     age                 -0.02705  0.01197   -2.26063  0.02383     I(age^2)             0.00102  0.00066    1.56271  0.11819     genderMale:age      -0.02109  0.01426   -1.47927  0.13913     genderMale:I(age^2)  0.00230  0.00079    2.89264  0.00384      So now from the summary, we can get the estimated values for each variable, however we have no information about the distribution of the errors, so how could I simulate from this specific model. The behind reason for why I want to do it is that I want to test whether the results from the quantile regression is reliable at the extremes of the data range,


I have three diagnostics tests. Each test indicates classifies an individual as either 0=ill or 1=non-sick. I assume these test are not independent. **How can multiple binary diagnostic tests be best integrated to form an overall diagnosis?**


title: Solving a practical machine learning problem

text: I am currently doing my Phd in computational biology at Stanford. I get the data I need to answer the questions I am interested in. The data sets are sometimes "large" and these large problems take longer time periods to solve (a couple of days sometimes). That being said I was wondering how machine learning on extremely massive data sets works? Suppose google wants to solve $Ax = b$ where $A$ has 10 billion rows, finding any gradients seems prohibitive. If google actually ran these simulations for as long as it takes (my equivalent of a couple of days), the solution maybe worthless before it arrives. This problem will be accentuated while training neural networks or implementing more complicated methods. What are practical solutions to this problem? I have seen statements like "We pick representative samples...". This is an absurd statement in my opinion because when p >> n, nothing is representative since the systems are under-determined. Any help on what 'representative' in these cases will also help.


title: How to prepare variables with mild skew for multiple regression?

text: I am doing some univariate analysis on a variable before doing regression. I think it is very skewed.   Three histograms are of (1) the original variable; (2) log10 transformation, and (3) inverse of the data to the power of 4.5. **How should the variable be prepared for regression?** ![enter image description here](http://i.stack.imgur.com/vIURC.jpg)


title: ANCOVA: ordinal covariate on SPSS

text: I'm trying to run an ANCOVA on SPSS with 3 variables: Gender, score 1, score 2. I want to control for the possibility that gender has an effect the scores, in order to isolate just the relationship between score 1 and score 2. However, I cannot enter gender as a covariate on SPSS as, by definition, it assumes covariates must be scalars. So, how can I run my ANCOVA and make sure that gender isn't affecting my 2 other variables?


title: Question about the error term in a simple linear regression

text: Suppose we have a linear regression model $Y_{it} = \beta_0 + \beta_1 X_{it} + \epsilon_{it}$, many times in literature it is assumed that $\epsilon_{it} \sim N(0,\sigma^2).$ This assumptions makes sense if we have a large data set due to the central limit theorem. My questions is that in certain situations I feel the error term being normally distributed is the wrong assumption. Suppose $Y_{it}$ is a bounded variable, such as age of a person, or a exam score of a student. Then if $\epsilon_{it} \sim N(0,\sigma^2)$ in this situations where $Y_{it}$ is bounded, is it not possible for the error term to be such that it forces $Y_{it}$ out of its bounds? For example suppose $Y_{it}$ represents a persons age, if the error term is normally distributed, then a random event could occur so it is possible for a person to live say a 1000 years? Hence, how do we fix this issue with the error term when our dependent variable on the left side on the linear equation is bounded. We could choose another bounded distribution for the error term, such as the uniform distribution over the bounds of $Y_{it}$. However this would not be realistic since it would imply all events in the error term are equally likely to occur. I am interested to here people thoughts about this problem. Edit: From reading all the great answers and comments below, here is what I have to say. Would it be practical to impose a bounded domain distribution on $\epsilon_{it}?$ For example the triangle density over a particular domain that $Y_{it}$ is in. Would imposing these types of distribution which have a bounded domain and resemble the normal distribution have any disadvantages?


title: How to measure whether the YES/NO (Bernoulli trials) answers have converged about their estimated mean for a specified error bound?

text: I have a questionnaire with questions that require a Yes or a No answer (1 or 0 input) from participants. From the set of answers I want to use the mean, $\mu$, estimated from the data for the likelihood (proportion) of yes answers. I would like there to be a confidence interval about this mean. I want to continue to request answers from new participants until a percentage, eg. 95%, certainty that the mean lies within this confidence interval. I assume that the underlying distribution is not normal because the domain is $[0,1]$. What are some possible methods to use to measure whether the confidence interval of the mean from the data is less that a certain chosen value?


title: Compare datasets from two machines

text: I have two data sets from two machines about their cpu, mem for n intervals of time. Machine 1 = [(cpu1, mem1), (cpu2, mem2), (cpu3, mem3), (cpu4, mem4), .... (cpuN, memN)] Machine 2 = [(cpu1, mem1), (cpu2, mem2), (cpu3, mem3), (cpu4, mem4), .... (cpuN, memN)] I need to find out how similar is the usage pattern of the (cpu, mem) in two different machines. In order to conclude how similar are the usage purpose of two machines Any statistical technique or Ml techniques thats appropriate to my problem. Please suggest. Thank you for your time.


title: Difference-in-differences with no pre-treatment?

text: The typical difference-in-differences estimator (as fixed effects) fits a model of the form $$ y_{it} = \alpha_i + \delta T_{it} + X_{it}'\beta + \epsilon_{it} $$ where $T$ is some treatment that happens to $i$ at time $t$. The coefficient $\delta$ is identified from the jump between time periods when T goes from zero to one, essentially using as counterfactuals those that didn't get treated during that period, after controlling for unobservables that don't vary in time. Normally the (panel) dataset starts with everyone un-treated, and ends with some remaining untreated while others get treated. Alternatively, if everyone (eventually) gets treated, you can still include post-treatment data to improve statistical precision -- the $\delta$ is still identified from the time periods where some got treated and others didn't. My question: is it legit to fit a model where one group starts treated, the other group starts untreated, and then the untreated group gets treated? This is basically the mirror image of a situation in which one group stayed untreated and one group got treated -- we still have heterogeneity in some time periods. Mathematically it seems identical -- the standard error components motivations seems to still apply. Am I missing something?


title: Bayesian stats: trick to accept the null?

text: There's a lot to be said and read about this, but I haven't found a clear answer to this question: Bayesian statistics are said to 'penalize' vague hypotheses with weak priors, by giving more support for the null hypothesis.   Say the theory I'm interested in proving actually predicts that the null is true. I could cheat, and set up my model with a weak prior. This would unfairly bias my empirical evidence towards the null. How can this be prevented? Thanks! EDIT: I realized I was actually referring to Lindley's paradox, where as far as I understand) a precise null and an uninformative prior might bias towards the null, whereas frequentist statistics would reject the null. http://www.laeuferpaar.de/Papers/LindleyPSA.pdf


title: Insignificant block effect in two-way block design

text: I used a two-way ANOVA to analyze data collected in a randomized complete block design: Y~Treatment+Block, where blocks represented spatial units Whereas Treatment had a significant effect on Y, block effects were insignificant. Does this simply mean that all blocks were the same with respect to the response variable? If that is so, would you it be a good idea to use an ANCOVA to see how parameters collected to describe spatial units affect Y? Thank you for your answers!


title: column proportions comparisons for one variable

text: I have a one variable with three categories like marriage/divorced/widowed. The counts are: 40, 23, 7. I use SPSS in this case. I know how to perform chi_squared test for this, and it shows that the counts are not equal (this is my null hypothesis). And I am to show that the first category shows up more often than other two. So the question is, how can I perform multiple comparisons for proportions? I mean something similar to column proportions comparisons (spss requires two variables).


title: Standard error bars overlap but significance - estimated marginal means versus observed means

text: I'm running a Mixed effects model ANOVA with two fixed factors (condition, repetition) and one random factor (subject). Subsequently, a Tukey multiple comparisons test is performed. Now I'd like to plot the means and standard errors (SEMs) of the single conditions in a single error bar plot, and report the p values between the conditions. The problem: while in the Tukey test, I got significant differences and non- overlapping SEMs between certain means, for my plotted real/observed data the SEM bars overlap. This is now counterintuitive, since commonly you would assume that in the case of overlapping, the means are not significantly different. My question is:   * is the difference between estimated marginal means and observed means due to having a random factor in my model, or what is the reason for the discrepancy?   * how would you report the data? Would you still plot observed data with the p values and state that the p values are derived from the estimated model? Or would you plot estimated means and standard errors? Thank you! EDIT: I'm adding the multiple comparisons result for a sample case as well as the observed means and standard error plot in case this helps. ![Observed means and standard errors](http://i.stack.imgur.com/iAy0X.jpg) ![Multiple comparisons result](http://i.stack.imgur.com/hWuCy.jpg)


title: Log-Cauchy Random Number Generation

text: I need to draw random numbers from a log-cauchy distribution which has density: $$f(x;\mu,\sigma)=\frac{1}{x\pi\sigma\left[1+\left(\frac{ln(x)-\mu}{\sigma}\right)^2\right]}.$$ Can anyone help me out or point me to a book/paper that could show me how?


title: Goodness-of-fit value from orthogonal distance regression

text: I have two sets of data, $x$ and $y$, with poissonian distributions. I want to check if the relation between $x$ and $y$ is a proportionality $y = ax + b$, so I used some algorithms to do Bivariate Correlated Errors and intrinsic Scatter (BCES). It is the first time I use it. These algorithms return the parameters values, their errors, and a $2\times2$ covariance matrix. I am trying to read tons of numbers and definitions, often ambiguous and/or conflicting, and anyway not easy-to-understand. My goal is to understand how good is the fit, that is if the line fits the data points. Something like the $\chi^2$ goodness-of-fit. Why the algorithms return the covariance matrix (instead of a goodness-of-fit value)? Is it possible to infer one from the other?


title: Replication reliability analysis

text: does anyone know about replication reliability analysis? I'm conducting a research which is similar to the study below (measure inter- and intra-subject reliability): http://link.springer.com/article/10.1007/PL00009575 May I know how can I analyze the data using replication reliability using SPSS? I've came across intra-class correlation (ICC) but I'm unsure if this is the correct statistical analysis. I've posted a question about ICC here: Inter- and intra- reliability (intraclass correlation) Thank you. Added: I would like to measure the inter- and intra-subject reliability in preparing three different samples (i.e. mild, moderate & extreme). Each subject is asked to prepare 9 samples (i.e 3 for mildly-thick, 3 for moderately-thick, and 3 for extremely-thick). Here's the data I get: ![enter image description here](http://i.stack.imgur.com/xT41C.png) How should I analyze the data in order to identify the inter- and intra-subject reliability. Previous study used replication reliability analysis but I can't find information about it.


title: Q learning in a stochastic environment

text: Most examples I have seen about Q learning, are performed in a deterministic world. For example, in the traditional grid world, the agent can finally do the path searching by exploring and exploiting the environment with a reward function without knowing the transition probability function. $$Q(s,a) = Q(s,a) + a*[ Reward + discount * Max Q(s',a') - Q(s,a)] $$ Now suppose the grid is a stochastic environment, an agent can move up/left/right with 1/3 probability. How can I program the Q learning, does that mean that in calculating the $Max Q(s',a')$,   $$Max Q(s',a') = Max [ P(up)*Q(s',up) , P(left) *Q(s',down) , P(right) * Q(s, right)]?$$


title: Can a nuisance multi-class classifier do better than binary classifier?

text: This is rather a theoretical question in order to save the trouble in trying to do empirical testing and is part of a bet, so I hope I am right... Say there are M classes in the data BUT you want to classify JUST between to subsets of these classes: M1 and M2 (M1+M2 = M). For example all M1 are different types of fraud and all M2 are different types of genuine users. Notice that the type of fraud/user is nuisance information and only the label fraud/genuine is of importance. One approach would be to use a multi-class classifier and then see if the estimated class is within M1 or M2. An alternative would be to use a binary classifier, disregarding the multiple labels and just using label = 1 for M1 and label = 2 for M2. Which classifier will work better in the general case? If the answer depends on data distribution, please explain. Thanks, Hanan p.s. my intuition says that binary will work better: the hypothesis space is smaller, so generalization error is smaller too.


title: Why does path analysis require variables to be measured on an interval scale?

text: According to this document: "In order to “correctly” use path analysis, Kerlinger and Fedhazur (1973) incidate that the following four assumptions must be met: (a) relations among variables are linear, additive, and causal (curvilinear and multiplicative relationships and inter— actions are excluded); (b) residuals (to be introduced later) are not correlated; (c) a one-way causal flow exists; and (d) variables are measured on an interval scale." What is the reason for this requirement of interval scale?


title: Testing a 2D point cloud for banana shape

text: I have many point clouds of a small size (say <200 points) in 2D. Some of them are isotropic and can be modelled as a single point. Others are elongated and curved so that they can reasonably be fit to the apex of a parabola. The dispersion seems to be of the same order for all clouds and there can be a fraction of outliers (say 10%). Here are four samples, of diminishing size. The bottom right is an isotropic case. ![enter image description here](http://i.stack.imgur.com/0ckvy.png) ![enter image description here](http://i.stack.imgur.com/ilkCN.png) ![enter image description here](http://i.stack.imgur.com/v0sYg.png) ![enter image description here](http://i.stack.imgur.com/N8Euq.png) I am trying to distinguish the circular clouds from the parabolic ones. I can probably use the regression coefficient, but I am looking for alternatives, with two purposes: 1) have a low computational cost (as the test has to be repeated hundred thousand times), and 2) go as far as possible as regards the classification of banana vs. "cherry" for small clouds or for larger dispersion, i.e. distinguishing cases like 3 and 4. The method can require the knowledge of the dispersion, but it is better if not.


title: sample size to estimate the best k in kNN

text: I would like to use a kNN classifier. My data set is quite small and include 2 classes having about 200 samples each one. I need to estimate k by using a cross-validation approach. How many samples could I use to find the best k?


title: information theoretic model selection

text: I have some questions on the information theoretic approach to linear modelling:   1. Why not just use residuals instead of likelihood?   2. Can same variables with different subsets of data be used to build a 'best model'?   3. Can same variables with different error distributions (for example poisson, negative binomial) be used to build a 'best model'?   4. Can same variables with different transformations (for example square root, log) be used to build a 'best model'?   5. Can interactions be included in candidate models that aren't included in global model?


title: Calculating possible combinations with restrictions

text: I have five additives that can be mixed into a chemical. Each one can be mixed in at a discrete percentage from 0.02% to 0.22% of the chemical (i.e. 0.02%, 0.03%, 0.04%,...0.22%). Each additive must be present, so mixed in at least 0.02%. The restriction is that the sum of all percentages cannot exceed 0.3%, but the sum can be less than that, so long as each additive is present. Without the restriction the answer is easy. But I can't figure out how to eliminate the number of combinations that would exceed the maximum allowed. If one additive is mixed at 0.22%, then the other four must be 0.02% each and only that. Thanks!


title: Contrastive Divergence: making hidden states binary

text: When learning a Restricted Bolzmann Machine with Contrastive Divergence $\mathbf{CD}_n$ for $n>1$ should we make the hidden states binary in the updates that are in between? In Hinton's - 'A Practical Guide to Training Restricted Bolzmann Machines' two contradictory reccomendations are mentioned: > 3.1 "It is very important to make these hidden states binary, rather than > using the probabilities themselves. [...] For the last update of the hidden > units, [...] use the probability itself" > > 3.4 "When [the hidden units] are being driven by reconstructions, always use > probabilities without sampling." The first quote seems to say the updates in the middle should do sampling (making binary), while the second seems to say they shouldn't! Which one is correct?


title: Test to show consistency in behaviour over an observation period

text: I am very new to the field of using Stats, and had some doubts regarding what test to use to analyse my data. So I have a set of individuals. I have observed these individuals for a particular activity over say 3-4 days. Now I want to say that the activities of these individuals remain consistent over days. I have used the Kruskal-Wallis test to show that the values on the different days are not significantly different. But I have to do something more, which is to show that the activity of each individual remains more or less consistent over the observation period. Does anyone have an idea about what statistical test I might use for this?


title: Calculating False Acceptance Rate for a Gaussian Distribution of scores

text: I have a biometric system that outputs a distribution of scores that resembles a Gaussian distribution (similar to the example graph in the following link: LINK). My point of confusion is how I calculate the False Acceptance Rate. How does threshold factor into the whole problem?


title: More than one function of a single instrumented endogenous variable -- is the model still underidentified?

text: **The typical instrumental variable setup** seeks a consistent estimate of $\beta$ from $$ y = \alpha + \beta x + \epsilon $$ where $cor(x,\epsilon) \neq 0$, in the univariate case, without loss of generality, by instrumenting with some variable $z$ that is uncorrelated with $y$ conditional on $x$ (the "exclusion restriction"). $\hat\beta_{IV} = (Z'X)^{-1}Z'y$. It can also be motivated as two-stage least squares, etc. In the multivariate case, one generally needs at least as many instruments as endogenous variables. **I've got a somewhat exotic case** , involving nodes in a network. $y_i$ is an outcome for node $i$, and $x_i$ is an endogenous predictor of $y_i$. I've got a single $z_i$ that is randomly assigned, and which predicts $x_i$ but not $y_i$, except through its effect on $x_i$. It gets tricky when I consider that $y$ is also predicted by some $f(x_{-i})$ (i.e.: values of $x$ at other nodes). The first stage model is something like $$ x_i = f(z_i) + \eta_i $$ where $cor(z,\eta) = 0$ The second stage, represented as two-stage least squares, is something like $$ y_i = \alpha + \beta_1\hat x_i + \beta_2 \hat u_i + \epsilon_i $$ where $u_i = f(x_{-i},D_{-i})$ -- basically some function of the $x$'s of nearby $i$'s, and their distances $D$. $\hat u$ is therefore $f(\hat x_{-i},D_{-i})$ So, I've got more endogenous variables than instruments, but my endogenous variables are all derived from a single endogenous variable, for which I have one instrument. **My question** : is my second stage underidentified? Or just identified? And why?


title: Is it possible to build a curve with a given skewness and kurtosis?

text: I have data in the following form:               Day  Kurtosis    Skewness   Total Inflows     163     .3          .5        435670      I'm attempting to do some modeling of an inflow regime between a river and a reservoir. To do this, I need a simulated inflow regime based on observed data. My original plan was to build a distribution based on a given kurtosis and skewness since I'm able to randomize those easily. The idea was that I could generate a curve from those two data points and the area under the curve to give me valid curves similar to the observed data. The resultant curve should look something like this (although there is significant variation year to year): ![Average inflows by day, with the green area indicating the portion I would like to simulate.](http://i.stack.imgur.com/ISTuk.png) I can add the daily variation after the curve is created without issue (I think). Unfortunately, I'm very new to the math and stats that are in play here, and I'm now unsure if it is possible to do anything with this data to give me fabricated curves. Does anyone know of a method for creating curves in this way? Alternatively, are there other methods that I might be able to employ to the same effect given different data? I also have data for each (observed) year in the following form:               Date             Inflows     12-Oct- 2002     23000      with a similar entry for each day. Thanks in advance for any insights!


title: difference in results of mixed effects model in nlme and lme4 package

text: Can anyone please tell me why the results of random slope model is different for the same dataset when I use lme and lmer. I first fitted a random intercept model as follows using both lme as well lmer               mdl1<-lmer(yld.res ~ rain + (1|state),data=data)   #random intercept model using lmer     mdl2<-lme(yld.res ~ rain,random= ~1|state,data=data) #random intercept model using lme          coef(mdl1)      $state      (Intercept)      rain     a   -336.4329 0.2711834     b   -294.2122 0.2711834     c   -256.1548 0.2711834     d   -263.4723 0.2711834     e   -217.1181 0.2711834     f   -239.2984 0.2711834          coef(mdl2)      (Intercept)      rain     a   -336.4333 0.2711836     b   -294.2125 0.2711836     c   -256.1550 0.2711836     d   -263.4726 0.2711836     e   -217.1183 0.2711836     f   -239.2986 0.2711836      As you can see, both of these yield the same results But when I try to fit a random slope model, both give different results:               mdl3<-lmer(yld.res ~ rain + (rain|state),data=data)     mdl4<-lme(yld.res ~ rain,random= ~rain|state,data=data)     coef(mdl3)     $state      (Intercept)       rain     a   -124.4119 0.09613782     b   -126.0181 0.11115529     c   -590.5186 0.65422357     d   -334.9604 0.35443209     e   -477.2628 0.61681345     f   -556.7407 0.65785116          coef(mdl4)          (Intercept)        rain     a   -16.09254  0.01476100     b    12.14178 -0.01015262     c  -761.00684  0.83513050     d  -327.16451  0.35018331     e  -451.16396  0.58277976     f  -632.17825  0.74185372      On top of that, when I run mdl3 (random slope using lmer) it gives the following message:               Warning messages:     1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :      Model failed to converge with max|grad| = 1.30318 (tol = 0.002, component 1)     2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :      Model is nearly unidentifiable: very large eigenvalue      - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio      - Rescale variables?      I am not sure why is this happening. Which random slope model should I use? The one with lme or lmer?


title: Why are 'random' and 'repeated' in mixed models in SAS both producing the same result?

text: Why does SAS random and repeated both produce the same result? Can someone explain this in detail? For example:               proc mixed data=test;     class variable1 ..... variableN;     model outcome=variable1+...+variableN;     random intercept/ subject=cluster type=cs;     run;          proc mixed data=test;     class variable1 ..... variableN;     model outcome=variable1+...+variableN;     repeated/ subject=cluster type=cs;     run;      Why do both program produce the same result?


title: robust standard errors newey-west in probit model

text: In R there are a way to calculate robust standar error newey when we adjust a linear model, using The function `NeweyWest()` in the `sandwich` package, but if you have a probit model how can do this ? there exist any function o you have to create a program. thanks


title: Module Updater Rejection

text: I was trying to build a Poisson model in OpenBUGS with the following codes:               model{     for(i in 1:N){          Y[i]~dpois(mu[i])     log(mu[i]) <- b[1]+b[2]*SU_WI[i]+b[3]*VMT[i]     }     for(j in 1:3){     b[j]~dnorm(0,0.001)     }     }      When I run these codes in OpenBUGS, the log file shows the following: > Sorry something went wrong in procedure Updater.Mode in module > UpdaterRejection Could you help me to solve this problem? I also need to build a multivariate crash model (3 dependent variables) in OpenBUGS with 5 independent variables where crashes are discrete events (Poisson-distributed). Can you provide me with some sample code to build such a model?


title: How to implement data I have to svmtrain() function in MATLAB?

text: I have to write a script using MATLAB which will classify my data. My data consists of 1051 web pages (rows) and 11000+ words (columns). The first 230 rows are about computer science course (to be labeled with +1) and remaining 821 are not (to be labeled with -1). I am going to label few part of these rows by myself. Then SVM will label the remaining unlabeled rows. I have found that I could solve my problem using MATLAB's `svmtrain()` and `svmclassify()` methods. First I need to create `SVMStruct`.               SVMStruct = svmtrain(Training,Group)      Then I need to use               Group = svmclassify(SVMStruct,Sample)      But the point that I do not know what `Training` and `Group` are. For `Group` Mathworks says: > Grouping variable, which can be a categorical, numeric, or logical vector, a > cell vector of strings, or a character matrix with each row representing a > class label. Each element of Group specifies the group of the corresponding > row of Training. Group should divide Training into two groups. Group has the > same number of elements as there are rows in Training. svmtrain treats each > NaN, empty string, or 'undefined' in Group as a missing value, and ignores > the corresponding row of Training. And for `Training` it is said that: > Matrix of training data, where each row corresponds to an observation or > replicate, and each column corresponds to a feature or variable. svmtrain > treats NaNs or empty strings in Training as missing values and ignores the > corresponding rows of Group. I want to know how I can adopt my data to `Training` and `Group`? I need (at least) a little code sample.


title: Obtaining SAS experience

text: I have been looking at analyst job postings and one of the most common requirement is experience of SAS.   * Unless your organisation currently uses SAS, how can you train as a SAS user?   * What programming language would be equivalent to SAS that employers might be happy to accept?


title: How do you create variables reflecting the lead and lag impact of holidays / calendar effects in a time-series analysis?

text: I am working on a time-series project in which I am forecasting the daily activity of something (let's call it 'Y') based on three years of historical data. I know that Y is affected by calendar effects such as Holidays (e.g., it will go down on Christmas). When I analyze the data, I want to account for the calendar effects. It is clear-cut how to create dummy variables for Holidays (e.g., 1 for Xmas and 0 for the remainder of days). It is less clear, however, how to create variables that account for the lead and lag effects of different calendar days. This is particularly important because a few of my days of interest (days that are not normal holidays but of interest to my study) are likely to have significant lead AND/OR lag effects. For example, activity should start increasing over the course of 30 days (30 is a simple estimate based on visual examination of the data) up to a particular calendar day of interest and for another calendar day of interest, activity should go up on that day and then remain up for a few weeks after, then slowly coming back down. In this context, I have two questions. One, is there a good statistical approach that I can use to determine what the lag and lead impact is of certain calendar days (i.e., how many days around the calendar day of interest appear to be impacted by that calendar day)? Any particularly open-source package in either python or R? Second, what should variables that incorporate lead and lag effects look like? For example, would a variable with a lead effect that builds up then drops off look like this: [0, 0, 0, 0, 1, 2, 3, 4 5, 6, 0, 0, 0]? In my first go-around with these data, I used a GBM with variables encoding holiday/calendar effects. For the calendar days that have substantial lead and lag impacts (impacts identified through visual examination and domain knowledge), I just created a variable that has days until that calendar day (e.g., [4, 3, 2, 1, 0, 365, 364, 363, etc.]. This model works fairly well but I think that I can improve it by improving how my variables represent calendar effects. Thanks a bunch for any responses.


title: Sum of Dependent Poisson Random Variables

text: I am working with the distribution of the sum of two dependent random variables. In my problem, there are two unobserved events, X and Y, where X precedes Y and Y is a function of the outcome of X, but I only observe the sum of the two outcomes, i.e: $$Z=X+Y|X$$ $$Pr(Z=x)=\sum_{z=0}^xPr(X=z)*Pr(Y=x-z|X=z)$$ It looks like a standard convolution distribution, except that Y and X are not independent. In my particular case, X and Y have a Poisson distribution, where the mean of Y is an increasing function of the value of the outcome of X. I am unaware of a closed form for the representation of the PMF of this distribution, and so need to derive it directly through the above sum. For my particular problem, I need to compute the PMF thousands of times, and am looking for ways to reduce the computational burden of doing so. As such, I was hoping someone would either 1) Know if there is a closed form representation for the PMF 2) Know of any approximations to the PMF Any help would be greatly appreciated. Thanks!


title: How to deal with unbalanced data

text: I'm doing data analysis with a dataset of 11795 data points (with 88 features). 85% (9973 points) of these data points correspond to data points belonging to class 1, 5% (589 points) belong to class 2 and 10% (1233 points) belong to class 3. I'm trying to build a model from this data for predicting the class of new data points. I started to wonder if I build my model using this dataset, does it favour the class 1 data points? Would it be difficult for the model to detect the low frequency classes? Generally how does one tackle unbalanced data sets such as the one I have? Thank you for any advices =) P.S. I'm using k-nearest neighbor and regularized linear regression methods.


title: Can I use averages of a likert scale as a dependent variable in ordinal regression?

text: I found similar threads on the forum, but none seems to apply to my case, so I decided to post the question. My response variable is an ordered 1-4 likert scale. I asked three subjects to rate the linguistic quality of sentences pre and post a linguistic revision process carried out on the sentences. I decided that the best way forward is analysing the post scores, having pre scores as a predictor in ordinal regression. Since I have three versions of the scale (3 subjects), if I use mixed-effects analysis with a random effect of subject, I'd end up having the same sentence with different pre scores in the data (from different subjects), which does not seem coherent. Would it be OK to use averages (both for the dependent variable and for the pre score)? I would still treat the dependent variable as ordinal (now with 9 levels), instead of being 1-4, it'd have intermediate categories, say 1, 1.3, 1.6, 2...


title: How to test if two populations have the same demographics?

text: I have 2 population samples for which I want to assess whether they have comparable demographics. What is the right way to test this, that includes multivariate combinations ? I.e I would like to check that both datasets have:   * the same proportion of men   * the same proportion of age [0-18]   * the same proportion of men aged [0-18] etc. Here is an example of the data. Ideally, I would like to do the test in R. ![count data of population demographics](http://i.stack.imgur.com/xqawP.jpg)


title: Relative importance of variables?

text: The following output is given: ![t1](http://i.stack.imgur.com/cjWmV.png) The task is to state, which variables, among thoses that are statistically significnat at 0.05, have the greatest and least relative importance on the fitted model? **First question** : I think I have to look at the absolute value of the standardized estimates to get an approximate ranking of the relative importance of the input variables on the fitted model, right? **Second question** : How can I calculate these values? I thought about dividing the estimate by the standard error, e.g. $-0.9699/0.0385$ but this is not equal to -0.2074?


title: Is it ok to here test statistical significance using permutation test?

text: I'm measuring electric signals from humans, who I have put into two different categories A and B. I would like know whether there is some difference between the two groups. My total sample size is $N = 14$ (seven subjects in both groups). Now, I do something to the signals and in the end get one number for each subject. How to test whether the two groups are statistically different? I was advised to do a permutation test (how?) because I don't know anything about the distributions of the numbers I'm comparing. Makes sense?


title: Mixed effects models issue and correlation question. Getting desperate

text: I have been trying to get answers for a while, and have looked up every help section and tried other stats sites. I have a thesis that I need to be completing shortly, but one set of data has been completely out of my league due to a small sample size. Basically, I have four broods (birds with chicks - precocial so they can peck at hatch e.g. chickens). I measured the proportion of time each brood fed every other day (approximately) until they flew. I had three broods in the ocean front and one brood in the mudflat habitats. Each data point within a brood is non-independent, so the only way to compare the mf brood against the OF broods in terms of time spent feeding (and other behaviours) would be to use the mean and use a Kruskal-Wallis or a one-sample Wilcoxon (mu=mean for MF brood). Of course the sample is too small to get anywhere near significance even though there is a clear pattern of MF feeding more and being less disturbed when I look at all the data points. It was suggested that I could use a mixed effects model so that each data point for a brood could be used and each brood would be considered a block. The person typed out some ideas for models but no one could tell me how to do it step-by-step and I don't really understand them. I found a source which would let me do it in Minitab (https://onlinecourses.science.psu.edu/stat502/node/72) but when I tried it it told me that it is unbalanced and I think it is because the first block is in one 'treatment' and the other three in the other 'treatment' but there is no overlapping of either. Another issue is that now I may have figured what the models do, I think they just try to explain the reasons why things are as they are, but I would still like to be able to say that MF feed significantly more than OF broods and that MF brood is disturbed less often. The model would only try to explain the reasons why rather than if they are significantly different (which is all I am really trying to figure out). If there is an easier way to compare two groups that have independent and non- independent data then I would greatly appreciate it. I actually thought of using a point biserial correlation using the means with MF coded as 0 and OF coded as 1. It works but is it valid? (http://faculty.vassar.edu/lowry/pbcorr.html) Is there anyone out there who could possibly help me with this? I am reaching my deadline and know that if I could just do this part then the rest should be okay. I am also a foreign student and my professor has said he can't help me long distance. I also have flush rates vs number of days to hatch. I wanted to see if birds were less likely to flush off their nest closer to hatch so used correlation. I have data for 7 nests and plotted each data point (n=77). The issue again is that not all points are independent but I need to show the flush rate for nest age (and not all nests were measured at the same age). Should I do a correlation for each nest and then is there someway of combining them to give me a p-value, or should I leave it as it is? Thank you SO much. I am sorry if my questions were too long. Rachel


title: Normalization/comparison-when the observed cardinality is a random variable, over finite sets

text: Given, $p_{ij}=\frac{|A_i|+|A_j|}{|A_i\cup A_j|}$ for sets $A_i$,$A_j$ $\forall i\not=j ,\forall i \in \left \\{ 1,2...n \right \\} $ and given the fact that $|A_i \cap A_j|>0$ for all $A_i,A_j$: it is clear that $p_{ij}$ > 1.   $p_{ij}\in\mathbb{R}^+, \forall i,j$ and $|A_i|\in\mathbb{Z}^+$ are random variables (observed cardinalities over finite sets are random variables). The $p_{ij}'s$ are calculated from an observed sample of $|A_i|$,$|A_j|$ and the cardinalities $|A_i \cap A_j|, |A_i \cup A_j|, \forall i,j$, are also known. Properties: Since, $|A_i|+|A_j| >|A_i|\cup |A_j|$ , $$\frac{|A_i|+|A_j|+\epsilon}{|A_i\cup A_j|+\epsilon}<\frac{|A_i|+|A_j|}{|A_i\cup A_j|};\forall \epsilon>0$$ i.e, for a fixed $|A_i|$ greater values of $|A_j|$ will always lead to smaller $p_{ij}'s$ and smaller values of $|A_j|$ will lead to greater $p_{ij}'s$ as for a fixed $|A_i|$, the numerator and the denominator increase or decrease by the same amount. Also,the $max(p_{\mathbb{.}j})$ for different columns $j$ can be different. Applied Statistical Context(a): The larger the $p_{ij}'s$ in a column $p$ compared to column $q$, makes us believe that the set $p$ is better than set $q$ purely in an applied scenario when I say 'better'. Applied Statistical Context(b): These are slightly similar to this applied example: I get 96/100 in a math test, and my brother gets 450/500. Who did better? There are ways to go about this if the observed proportions are considered binomial with sample sizes 100, 500 respectively. Here, our sample size in the denominator is fixed, unlike my example where the denominator is a random variable. Now in order to analyze the columns of the matrix $P$, with entries $p_{ij}$, a primary issue is that the range of the values in different columns are not the same. How can I transform/(also may be normalize) the columns such that they are all within a fixed range [l,u] and also the entries in different columns need to taken to be account, instead of scaling individual columns. This may require tweaking my definition of $p_{ij}$ or might require some transformations over the columns taken together -or both. Also, feel free to make distributional assumptions over $|A_i|$-if you would prefer that direction.


title: Regression by means of polyserial correlations in R

text: I have two continuous independent variables and a binary dependent variable, and I want to use regression by means of polyserial correlation, but I can't find any information on how to do this in R. I know the alternative is to use logistic regresssion, but the problem that I'm trying to solve necessitates that I use polyserial correlation. I'm particularly concerned with getting the residuals of the regression. Any help would be much appreciated. If you know of another software that can perform this regression, I would also like to hear about it.


title: Which of these points in this plot has the highest leverage and why?

text: I am studying the definition of leverage, and I understand it in terms of formulas. However, if I would have a plot like this for instance, how could I see which of these points has the highest leverage? Which one would it be in this plot for instance?![enter image description here](http://i.stack.imgur.com/CwBtb.png)


title: How to get odds ratio when combining multiple tests in logistic regression in spss?

text: we get odds ratio in logistic regression for predictive accuracy of a test. how do we get one odds ratio from combining more than 1 test in predicting outcome? how do we enter the data to get this in spss?


title: Fit nonlinear parameter

text: I'm attempting to fit this model: $P = C_0 + C_1*U^r$ Given known vectors of observations $P$ and $U$, I want to fit values for $C_0$, $C_1$ and $r$. How do I make this fit in R? or preferably GSL in C? Is it possible?


title: Generalized Least Squares vs Ordinary Least Squares under a special case

text: This question regards the problem of Generalized Least Squares. Vectors and matrices will be denoted in bold. **Premises**. Let $N,K$ be given integers, with $K \gg N > 1$. The transpose of matrix $\mathbf{A}$ will be denoted with $\mathbf{A}^T$. Suppose the following statistical model holds $$ (*) \quad \mathbf{y} = \mathbf{Hx + n}, \quad \mathbf{n} \sim \mathcal{N}_{K}(\mathbf{0}, \mathbf{C}) $$ where $\mathbf{y} \in \mathbb{R}^{K \times 1}$ are the observables, $\mathbf{H} \in \mathbb{R}^{K \times N}$ is a known full-rank matrix, $\mathbf{x} \in \mathbb{R}^{N \times 1}$ is a deterministic vector of unknown parameters (which we want to estimate) and finally $\mathbf{n} \in \mathbb{R}^{K \times 1}$ is a disturbance vector (noise) with a known (positive definite) covariance matrix $\mathbf{C} \in \mathbb{R}^{K \times K}$. The Maximum Likelihood (ML) estimate of $\mathbf{x}$, denoted with $\hat{\mathbf{x}}_{ML}$, is given by $$ (1) \quad \hat{\mathbf{x}}_{ML} = (\mathbf{H}^T \mathbf{C^{-1}} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} \mathbf{y} $$ and this is also the standard formula of Generalized Linear Least Squares (GLLS). Consider the standard formula of Ordinary Least Squares (OLS) for a linear model, i.e. $$ (2) \quad \hat{\mathbf{x}}_{OLS} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T \mathbf{y} $$ As a final note on notation, $\mathbf{I}_K$ is the $K \times K$ identity matrix and $\mathbf{O}$ is a matrix of all zeros (with appropriate dimensions). * * * Now, for the problem at hand, assume that $\mathbf{C}^{-1} = \mathbf{I}_K + \mathbf{X}$, where $\mathbf{X} \in \mathbb{R}^{K \times K}$ is a symmetric, invertible matrix (and, for the formalism to make sense, $\mathbf{X}$ is such that $\mathbf{I}_K + \mathbf{X}$ is invertible and the inverse is positive definite). In this case, it can be proven (using _matrix inversion lemma_ ) that $$ (3) \quad (\mathbf{H}^T \mathbf{C^{-1}} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T + \mathbf{Q} $$ where the expression for the matrix $\mathbf{Q} \in \mathbb{R}^{N \times K}$ can be found (I will omit it here). Now, finally, **Proposition 1**. _If $\mathbf{H}^T\mathbf{X} = \mathbf{O}_{N,K}$, then equation $(1)$ degenerates in equation $(2)$, i.e., there exists no difference between GLLS and OLS_. The proof is straigthforward and is valid even if $\mathbf{X}$ is singular. Now, my question is **_Question_** : _Can an equation similar to eq. $(3)$ (which "separates" an OLS-term from a second term) be written when $\mathbf{X}$ is a singular matrix?_ I am not interested in a closed-form of $\mathbf{Q}$ when $\mathbf{X}$ is singular. But I do am interested in understanding the concept beyond that expression: what is the actual role of $\mathbf{Q}$? In which space does it operate? I found this problem during a numerical implementation where both OLS and GLLS performed roughly the same (the actual model is $(*)$), and I cannot understand why OLS is not strictly sub-optimal. I found this slightly counter- intuitive, since you know a lot more in GLLS (you know $\mathbf{C}$ and make full use of it, why OLS does not), but this is somehow "useless" if some conditions are met. What are these conditions? As a final note, I am rather new to the world of Least Squares, since I generally work within a ML-framework (or MMSE in other cases) and never studied the deep aspects of GLLS vs OLS, since, in my case, they are just intermediate steps during the derivation of MLE for a given problem. If the question is, in your opinion, a bit too broad, or if there is something I am missing, could you please point me in the right direction by giving me references? Preferably well-known books written in standard notation.


title: logistic regression with nominal clustered responses

text: If we would like to study the admission rate with gender as explanatory variables in certain university, we generally cluster the result within each independent department. We would generally assume there is nonnegative correlation between gender and department. for the logistic regression model $logit(P(Y_i=1))=\mu_i+\beta x_i$, where $x_i$ is indicator variable for gender. There is no correlation between gender and department in the model above. How do we quantify the correlation into the model without introducing two many explanatory variables? If we stick to the logistic model above, Why we could assume the gender effect is the same among different department?


title: Why cov(AX)=A cov(X) A'

text: I cannot verify the following theorem. Maybe I am doing something wrong, but I don't know what?! Additionally, I'm not sure about the meaning of a constant matrix in the theorem. Theorem: cov(AX)=Acov(X)A' if A is a constant matrix. (A' means transpose of A, i.e., A^T) Example:               A (10x5):              0.8727    0.7145    0.9147    0.6370    0.9760         0.8241    0.6210    0.0735    0.7318    0.1231         0.8609    0.7410    0.9865    0.6090    0.7693         0.1135    0.9563    0.2369    0.5369    0.9876         0.7349    0.4313    0.9650    0.1657    0.4418         0.7872    0.6589    0.4490    0.5974    0.7669         0.4247    0.6411    0.9335    0.5238    0.9153         0.2126    0.6212    0.6423    0.3651    0.9752         0.5092    0.6179    0.3732    0.9569    0.1673         0.8009    0.4691    0.4420    0.6289    0.3132          X (5x5):              0.6635    0.6323    0.5883    0.9362    0.0998         0.8298    0.3788    0.4760    0.8170    0.3296         0.8553    0.6836    0.4701    0.5723    0.6803         0.0065    0.8085    0.1542    0.5780    0.7653         0.6389    0.9297    0.1418    0.0438    0.0886          AX (10x5):              2.5820    2.8701    1.5202    2.3352    1.5188         1.2083    1.5127    0.9452    1.7493    0.9079         2.5253    2.7069    1.5259    2.3616    1.5355         1.7059    1.9482    0.8562    1.3768    0.9861         1.9542    1.8325    1.1795    1.7078    1.0380         1.9469    2.2502    1.1887    1.9111    1.1263         2.2004    2.4239    1.2045    1.7985    1.3707         1.8312    2.0105    0.9173    1.3278    1.0287         1.2828    1.7403    0.9404    1.7555    1.2554         1.5029    1.7859    1.0436    1.7633    1.0443          cov(X) 5x5:              0.1190   -0.0345    0.0457    0.0139   -0.0482        -0.0345    0.0429   -0.0316   -0.0557    0.0015         0.0457   -0.0316    0.0419    0.0563   -0.0132         0.0139   -0.0557    0.0563    0.1175    0.0102        -0.0482    0.0015   -0.0132    0.0102    0.1009          cov(AX) (5x5):              0.2227    0.1984    0.0988    0.1011    0.0827         0.1984    0.1960    0.0916    0.1063    0.0874         0.0988    0.0916    0.0576    0.0745    0.0455         0.1011    0.1063    0.0745    0.1144    0.0618         0.0827    0.0874    0.0455    0.0618    0.0510          but A*cov(X)*A' (10x10):            0.2195     0.12119     0.20872    0.077985     0.15076     0.16413     0.17298     0.11964     0.14558     0.15641      0.12119     0.09362     0.12556  -0.0017266     0.10372     0.08966    0.081364    0.034694    0.099912     0.11002      0.20872     0.12556     0.20256    0.057506     0.15142     0.15574     0.15915     0.10193     0.14583      0.1577     0.077985  -0.0017266    0.057506     0.11155    0.015719    0.059279    0.089875    0.099947    0.026385    0.015011      0.15076     0.10372     0.15142    0.015719     0.12378     0.11292     0.10434    0.055612     0.10669     0.12511      0.16413     0.08966     0.15574    0.059279     0.11292     0.12295     0.12914    0.090107     0.10648     0.11599      0.17298    0.081364     0.15915    0.089875     0.10434     0.12914     0.14732     0.11382     0.11134     0.11072      0.11964    0.034694     0.10193    0.099947    0.055612    0.090107     0.11382     0.10483    0.061491    0.057124      0.14558    0.099912     0.14583    0.026385     0.10669     0.10648     0.11134    0.061491     0.12616     0.12169      0.15641     0.11002      0.1577    0.015011     0.12511     0.11599     0.11072    0.057124     0.12169     0.13273


title: Is it OK to take "Analysis of Variance in Experimental Design Models" Before a Regression Course?

text: I am a math major starting a MS program in statistics. I have taken an undergraduate sequence in stats (basics of distributions, estimation, hypothesis testing) but never a formal course in regression. There is a course called "Analysis of Variance in Experimental Design Models" which covers the following topics (Multivariate normal distribution; quadratic forms; linear models; simple random, randomized block, Latin squares, factorial, split-plot, balanced incomplete block designs; analysis of covariance; confounding; and multiple comparison tests.) which is more convenient for me to take before I would take the first course in regression. Is it OK to take this exptl design course _before_ I take a regression course or would it be best to take the class on linear regression first (so the explt design class would make more sense to me)? Thank you, JZ


title: Using machine learning to tell apart users from crawlers

text: Would it be possible to build anomaly detection algorithm to be able to tell apart real user from crawlers/spiders with high confidence? I have access to user tracking data where information about client session, IP address and weather client can process javascript. Based on that data I can calculate how many pages client visited and what is the entropy of time between requests. I also know which sessions converted which would be my non anomalous data I can use to train. To summarize I have 3 features I think I can use:               - number of requests     - *calculated entropy based on time spent on each page     - can run JS      Are there any other possible or better features that could be used? Is this problem solvable with machine learning? `*` crawlers usually don't spend too much time on the page so time between requests would be very short with similar intervals. That would give small entropy value. Real visitors would produce more random pattern.


title: Multiple testing: contingency table or logistic regression approach for survey?

text: I have a survey asking a person's gender (M/F), geographic region (West, East, Midwest, South) and a few other demographic variables along with two dozen questions in which they can disagree / slightly disagree / slightly agree / agree. The data was given to me and whoever worked on the data has refactored the responses to a simple disagree / agree. I am interested in whether people of different gender / region / other demographic variable have differences in responses. For example, Q1 could be something similar to 'You are stressed often' and the responses are either agree / disagree. I would want to know if gender is related to agreeing/disagreeing with Q1. The analysis suggested to me was to build a set of contingency tables for every combination of demographic vars and questions. For example: Gender vs. Q1, Region vs. Q1, Gender vs Q2 and so on. Then I would use Fisher's exact test or the chi-squared test of independence on each of the contingency tables. Example table: > >         Q1: Agree Disagree >     Male       10       30 >     Female     15       32 >   I have two dozen questions and at least 4 demographic variables. If I were to do all combinations I would have to run at least 96 chi-squared tests. Obviously this also brings in a multiple testing issue. As an alternative route, a colleague of mine suggested logistic regression with the following model: > Question Response ~ gender + region + other demographic vars + interactions. This would mean fitting slightly fewer models but will still have the same issue with multiple testing. I would still need to fit a logistic regression to each question. There's also the issue of interpretability: it's a lot easier to explain the results of a chi-square test than to explain what a logistic regression is. Which technique is preferable? Is there a better way to handle this data? In either case, how should I handle the multiple testing?


title: MOOC or other online class for advanced stats?

text: Are there any MOOC or other (free/cheap) online classes for advanced statistics? I'm thinking upper-level undergrad or introductory graduate level. Things like time series models, factor analysis, MLE, Bayesian statistics, etc.


title: Approach to post-hoc test for factors in a regressions (in R)

text: I have a problem with some analysis I need to do. I have a series of regressions. Some of the predictors of these regression are categorical with multiple levels. I performed regressions, both linear and logistic, choosing a baseline for these category according to various factors. The problem is that my colleagues asked not only for a confrontation of the factors to a baseline but also a pairwise confrontation. Like you it's used to do with a post-hoc test for ANOVA (they are pretty new to regressions and their benefits). How should I approach this? I thought of some solutions:   * Subsetting: That is subset the data to include two factors at time, and therefore repeating the regression once per every subset.   * Splitting: Splitting the category column in a column for every factor and put 0 and 1 as levels. This approach can furthermore be conducted in two ways:      * Putting all the new columns in the regression (minding that they are mutually exclusive).     * Putting one column at time, multiplying the regressions. Which approach would you suggest, minding statistical correctness and workload? Especially, what's the conceptual difference between the three methods? Thanks a lot!


title: Calculating the MacKinnon empirical distribution test to test mediation

text: I am interested in using the MacKinnon product of coefficients test to assess mediation. Is there a quick way to carry out these calculations and determine significance based on the obtained statistic?


title: Rank Order Analysis with Missing Values

text: I have several survey questions that allow a respondent to choose their top three choices from a set of potential answers--basically "What is the most important, 2nd most important, 3rd most important." This leaves me with a rank order question containing several missing values since there are more than just three options to choose from. I am aware of some of the traditional methods of analysis for rank order (i.e. Wilcoxon, Spearman, Friedman Test, etc.) but I believe that these all require every option to be ranked, not just a "top nth." Are there any suggestions on the best way to analyze this type of data?


title: What is the role/purpose of hidden variables in graphical models?

text: Is there a formal treatment of the role/power of latent/hidden variables in graphical models and other machine learning models (e.g., structural equation models)? For example, the Restricted Boltzman Machine is a graphical model which is a compact representation of a probability distribution of $x$, but it is modeled as a joint distribution with a hidden variable $h$ (i.e., we model $p(x,h)$ instead of $p(x)$). Why is $p(x,h)$ any more powerful than a suitable representation for $p(x)$?