title: Post-hoc tests for MANOVA

text: I am using a MANOVA test to compare nine different dependent variables (from neuropsychological and neuropsychiatric assessment) between three groups. The output shows a significant influence from GROUP on my variables ( _p_ < .001). Ofcourse, I am interested in how the three groups influence every dependent variable. I have studied Field's "Discovering Statistics Using IBM SPSS Statistics" chapter 16, and he states that the preferred post-hoc analysis is a discriminant analysis, because of the linear combination in which the dependent variables are related to group membership in a MANOVA. Discriminant analysis could account for this linear combination, so Field states. Otherwise, I read some literature, on basic statistical sites, where is stated that I can use multiple univariate ANOVA's with Bonferroni correction, and use a post-hoc on these univariate ANOVA's when they are significant. Which of these method's is better? That is, which will make my chance of a Type I or II error the least?


title: Performing logarithmic multiple regression with Excel?

text: Is it possible to perform logarithmic regression on multiple variables with Excel? If I just have a single independent variable than it's very easy to do this using the best-fit line option (it lets me switch from linear to logarithmic). But this feature does not work for multiple variable regression and the regression feature under the Data Analysis plugin only seems to support linear multiple regression. However, I have a table that has 3 columns containing 3 independent variables and 1 column with the corresponding dependent variable (outcome). I'm pretty sure there's a logarithmic relationship, but I'm not sure how to use Excel to get the coefficients. Thanks!


title: Cointegration-based feature selection

text: I have a basket of time series (stock prices). I want to find the N (fixed or not) time series that will best replicate the basket in the sense that combination of them will be best cointegrated with the basket.   Beside using the N series that have the best cointegration scores (ADF test) and regressing these variables against the basket. Do you know of other/better methodologies for doing that?


title: What will be minimum variance unbiased estimator?

text: Let $X_1, X_2, ..., X_n$ be a random sample from a distribution with p.d.f., $$f(x;\theta)=\theta^2xe^{-x\theta} ; 0<x<\infty, \theta>0$$ Obtain minimum variance unbiased estimator of $\theta$ and examine whether it is attained? MY WORK: Using MLE i have found the estimator for $\theta=\frac{2}{\bar{x}}$ Or as $$X\sim Gamma(2, \theta)$$So $E(X)=2\theta$ $E(\frac{x}{2})=\theta$ so can i take $\frac {x}{2}$ as unbiased estimator of $\theta$. I'm stuck and confused need some help. Thank u.


title: Comparing logistic regression coefficients across models?

text: I've developed a logit model to be applied to six different sets of cross- sectional data. What I'm trying to uncover is whether there are changes in the substantive effect of a given independent variable (IV) on the dependent variable (DV) controlling for other explanations at different times and across time. My questions are:   * How do I assess increased / decreased size in the association between the IV and DV?    * Can I simply look at the different magnitudes (sizes) of the coefficients across the models or do I need to go through some other process?    * If I need to do something else, what is it and can it be done/how do I do it in SPSS? Also, within a single model,   * Can I compare the relative size of independent variables based on unstandardised scores if all are coded 0-1 or do I need to convert them to standardised scores?    * Are there problems involved with standardised scores?


title: What quantifies a stable parameter?

text: I'm optimizing 5 parameters for an option pricing model. Now I want to asses whether these parameters are stable over time (i.e., a year). For this I create about 12 subsamples and estimate the optimal parameters for each subsample (I do this by minimizing a cost function, i.e., chi-squared, under a certain threshold). Now I'm not really sure what statistic quantifies those estimated parameters as stable. Is it a low standard deviation? What would be the best way to prove those parameters are stable over time? **EDIT:** The model which I'm using is neither differentiable nor is it continuous. Related question, without an answer: How to test stability/variability of statistic over time, using R


title: In search of a proper similarity function

text: I'm trying to find the most similar sample between a candidate and a bag of samples. Consider you have a knowledge corpus as follows:               corpus = data.frame(region = factor(), time = numeric(), observation = numeric())     corpus = rbind(corpus, data.frame(region = 'america', time = 1:20, observation = rnorm(20)))     corpus = rbind(corpus, data.frame(region = 'asia', time = 1:20, observation = rnorm(20)))     corpus = rbind(corpus, data.frame(region = 'europe', time = 1:20, observation = rnorm(20)))      Now I have a new sample and I want to choose from my knowledge base, which of the records is the most similar to my new sample. I used different metrics with not so much improvements between them:   * Euclidian    * Cosine   * Correlations But what I really want to see is that the matching is done based on the shape of the curve not in vectorial distance but I have no clues about similarity functions that would metric this type of space. *Edit: I want to match my new sample with the "closest" element in my knowledge base. But what I want is to have a match using the shape of the curve to compare the samples. If it has positive or negative derivatives in the same intervals, same amount of local maximums, that sort of things. I tried different approaches and compared the quality of them with a random assign using a KolmogorovSmirnov test. Fact is that my random assign is not showing to be worse than my super califragilistic detection method . Any clues? Thoughts?


title: maximizing frequency

text: I am trying to find the area with thickest density. I have newborn numbers by day for which I am trying to find the series of 7 consecutive days that would produce the highest frequency. So the start and end dates are not tied to any particular day of the week the period can start on any day but would end in 7 days from the start. I know I can iterate through the various starts but is there a more efficient mathematical solution?


title: Multidimensional scaling on distance or similarity matrix

text: Why doesn't the scatter plot change when I perform multidimensional scaling on distance or similarity matrix? This figure uses similarity matrix ![enter image description here](http://i.stack.imgur.com/ZKWYA.png) And this figure use distance matrix (sqrt(1-similarity)) ![enter image description here](http://i.stack.imgur.com/IVNCT.png)


title: Looking for name of statistical property (variance explained by regressor in addition to all other regressors)

text: There is a statistical property - and I only know the name in German and cannot find any proper translation - that is defined as: variance explained of a specific regressor x1 in addition to the variance explained by all other regressors. That is, the shared variance is _not_ included in the explained variance of this statistical property - it is the minimum estimate of its shared variance with the dependent variable. It's a classical concept in psychology. In German the term is "Nützlichkeit", which translates into "utility". But I'm afraid in English the name is something so different that my google-fu does not work.


title: Singular value decompostion covariance matrix (numerical recipes)

text: I am trying to implement singular value decomposition in C. I am using routines svdfit and svdvar from Numerical Recipes. The results from svdfit seem to be correct, but the results from svdvar are not. I am struggling with this problem for some time now, and I don't know how to solve it. For example, take the following data:               X = {1.47,1.50,1.52,1.55,1.57,1.60,1.63,1.65,1.68,1.70,1.73,1.75,1.78,1.80,1.83}     Y = {52.21,53.12,54.48,55.84,57.20,58.57,59.93,61.29,63.11,64.47,66.28,68.10,69.92,72.19,74.46}      For a 2nd order polynomial, the results should be:               a_0 = 128.8128 +- 16.3083     a_1 = –143.1620 +- 19.8332     a_2 = 61.9603 +- 6.0084      However, I have (the coefficients are obtained from svdfit and the standard errors are the diagonals of the covariance matrix obtained from svdvar):               a_0 = 128.8128 +- \sqrt{208104.0000}     a_1 = –143.1620 +- \sqrt{307785.6091}     a_2 = 61.9603 +- \sqrt{28247.8801}      Does anyone know what's the problem? Thanks so much for your help.


title: how do i interpret trace statistic of Johansen's cointegration test?

text: how do i interpret the results of the Johansen cointegration test if the result from my STATA output fails to accept both null hypotheses of r=0 and r<=1 ? for example: STATA Screenshot can i conclude whether there is a cointegrating relationship?? many thanks in advance!


title: Likelihood estimation using iid normal samples

text: Given an i.i.d. sample $X = (x_{1}, \dots, x_{n}) \sim N(\mu, 1)$. I have been asked to show that the likelihood of $\mu$ based on the whole sample is proportional to the likelihood based on $\bar{x}$ (the sample mean) alone. I'm a little confused by the phrasing of this question. I was hoping the perhaps somebody could elaborate on what this question is asking. Also, as an aside, given a random sample $X = (x_{1}, \dots, x_{n}) \sim N(\mu, \sigma^{2})$, how can one estimate the $\text{Pr(X > x)}$? Can the MLEs for $\mu$ and $\sigma^{2}$ be used in normal CDF to estimate $\text{Pr(X > x)}$? where $\hat{\mu} = \bar{x} = \frac{\sum_{i}x_{i}}{n}$ and $\hat{\sigma}^{2} = \frac{\sum_{i}(x_{i}-\hat{\mu})^{2}}{n}$ hence $\text{Pr}(X > x) = 1 - F(x)$ where $F(x)$ is the CDF evaluated at $x$.


title: How to Improve the Design of a Survey of Multiple Rare Variables Using Rare Variable Location Information

text: I am attempting to update the design of an annual survey which measures the percent-deficient of 26 different asset inventories. Most of these assets don't occur on every sampling unit and some of them are quite rare (occur on around 5% of sampling units). The asset deficiency rates can be estimated by either a simple proportion estimator, a domain proportion (for those that either exist or do not exist on a sampling unit and, if they exist, are either deficient or not), or a ratio estimator (for assets that have 0 or more units of inventory on a sampling unit and some or none of that inventory is deficient). The sample size that has been used historically is around 1,200 (~1%), which has been proportionally allocated into 72 strata. Most estimates from past samples have indicated useful accuracy (standard errors of 1% or less). However, for assets that are rare and have high deficiency rates the accuracy is much worse and not useful for their purpose (e.g. standard errors of 5%). Knowledge of which sampling units contain which assets is available (not necessarily the amount of that asset) and I would like to use this information to alter the sample inclusion probabilities within the same stratification scheme. Ideally, all assets would be estimated with a certain degree of accuracy. One way to do this might set inclusion probabilities such that the standard error would likely fall below some threshold for each asset (perhaps unique for different assets). If this is impossible for all 26 assets, then minimize the deviation above the threshold with more weight given to more important assets. My initial thought was to create a linear program to determine inclusion probability adjustment factors for a sampling unit based on which assets existed on the unit (e.g. if asset 1 was on a sampling unit increase the inclusion probability by some amount). A feasible solution would ensure that a sample's expected number of units which contain each asset would meet a requirement (e.g. at least 300 units containing asset 1 are selected). This requirement could be determined based on (1) the number of sample units containing the feature from the previous year and (2) the asset's standard error from the previous year. The LP objective might be to minimize the difference between the smallest and largest weights. Does this seem reasonable/possible? Are there better approaches? Are there references available?


title: Canonical parameter for Poisson Regression

text: I was looking at the Poisson regression. Here I find that the canonical parameter is taken to be the logarithm of the rate of the Poisson. I was wondering if there is any reason for considering this as the canonical parameter, rather than the rate of the poisson itself (with the restriction that the parameter is positive).


title: Will a quadratic and cubic latent growth curve model always have better fit than a linear model?

text: A conjecture that I have heard is this: _When you add a quadratic or a cubic term to a linear latent growth curve model, the fit will always improve_. Is this correct? Why/why not? As I have had the idea explained to me, adding additional terms will capture a additional portion of the variance and should therefore increase the amount of variance explained, therefore increasing the goodness of fit.


title: Estimated error variance $\sigma^2$ for MCMC estimation in a high-dimensional space

text: Let $f$ be a function such that: $$f~:~(x,~\theta)\in\mathbb{R}^{3}\times\mathbb{R}^{12} \rightarrow f(x,~\theta)\in\mathbb{R}^3$$ My observations $y$ are noisy values taken by the function $f(\cdot ,~\theta)$ for known values of $x$. I would like to estimate the conditional distribution of $\theta \in \mathbb{R}^{12}$ given the observations $y$. I assume the prior distribution of $\theta$ is uniform, and the likelihood of the observations given $\theta$ is Gaussian. I would like to have some insights about the estimation of $\theta$ with MCMC Metropolis-Hastings simulations. I use Gaussian proposals. To test the algorithm, I have created synthetic noiseless $y$ given a known $\theta$, and I try to use MCMC to find point-estimates really close to the true $\theta$. To make things harder, I start far from the true values, but even if I started "close" (in $\mathbb{R}^{12}$) to the solution, the convergence is not fast or obvious. I guess this is due to the quite complex non linear expression of $f$. 1) First of all, since $y\in\mathbb{R}^3$, would it be relevant to define separate Gaussian likelihood functions, one for each component of $y$, with potentially different error variances $\sigma_i^2$. Or is it better to just define one error variance $\sigma^2$ and jointly consider all the components of $y$ ? 2) Second, if I want to provide an estimate of $\sigma^2$ to the algorithm, I could just compute the empirical variance, which is roughly $\frac{1}{n-p}$ times the sum of square differences between $y$ and $f(x,~\theta_{init})$. However, I have noticed, by looking at the each simulated component of $\theta$ independently, that the chain does not mix well unless I give a much higher value to $\sigma^2$ (like omitting the factor $\frac{1}{n-p}$ where $n$ can be of the order of $10^5$ and $p=12$). Would this be due to the high dimensionality of the estimation problem? As I understand it, the higher $\sigma^2$, the higher the probability to accept a new move even if the move "looks bad", and thus the higher the acceptance rate. So, to put my second question in other words, should I be patient and drastically increase the number of simulations, or should I just force the chain to sample more space more quickly by increasing $\sigma^2$? I would love to get a result in $10^5$ simulations, but I am afraid I will need a lot more because of the dimensionality and non linearity.


title: R: forecast function accuracy for ARIMA models

text: I have a problem with the `forecast` function for ARIMA models in R. It calls `predict` that calls `KalmanForecast`. Ok...here's the deal. the mean one-step forecast of the Arima object produced by this call               forecast(Arima, h=1)$mean[[1]]      is often significantly different from the result of a manual forecast by conditional expected value (best linear predictor). For example a non seasonal Arima(1,1,1) without drift has of course the structure               y[t] = y[t-1] + AR1*(y[t-1] - y[t-2]) + MA1*epsilon[t-1] + epsilon[t]      so the one-step prediction is very straightforward               y[t] = y[t-1] + AR1*(y[t-1] - y[t-2]) + MA1*epsilon[t-1]      but this result is always different from the result of the forecast function call. Is it due to approximation errors in the Kalman recursion? Try yourself with this code, it only needs the `forecast` package                 x = arima.sim(n = 1000, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),sd=sqrt(0.1796))       t = length(x) + 1      Arimafit = Arima(x = x, order = c(1,1,1), seasonal = list(order = c(0,0,0), period =     1), include.mean = FALSE,include.drift = FALSE)      manualforecast = x[t-1] + coef(Arimafit)[["ar1"]]*(x[t-1] - x[t-2]) + coef(Arimafit)    [["ma1"]]*Arimafit$residuals[t-1]          autoforecast = forecast(Arimafit, h = 1)$mean[[1]]      `autoforecast` is always different from `manualforecast`, sometimes significantly.


title: Mixture distribution fitting for latent variable analysis

text: Are there any **analytic approaches** to using _mixture distribution fitting_ for _latent variable analysis_? I'm specifically interested in existing approaches to determining whether mixture components indicate a presence of _sub-populations_ OR _effect(s)_ from underlying _factor(s)_. References to relevant key _papers_ , _tutorials_ , _guidelines_ , _examples_ and `R` _packages_ are additionally welcome. P.S. This question implies the subsequent planned use of standard social sciences _quantitative multivariate methods_ , such as _factor analysis_ (`EFA` and `CFA`), in a target research study.


title: Is Student's t test a Wald test?

text: Is Student's t test a Wald test? I've read the description of Wald tests from Wasserman's _All of Statistics_. It seems to me that the Wald test includes t-tests. Is that correct? If not, what makes a t-test not a Wald test?


title: How to cluster time series?

text: I have a question about cluster analysis. There are 3000 companies, which have to be clustered according to their power usage over 5 years. Each company has values for every hour during 5 years. I would like to find out if some companies have the same pattern in usage power over the time period. The results should be used for daily prediction of power usage. If you have some ideas how to cluster time series in SPSS, please share with me.


title: Interpretation of null and alternative hypotheses

text: ![Question for interpretation](http://i.stack.imgur.com/nGx8c.png) The way I interpreted this question was that $H_0: \mu \ge 10, H_a: \mu < 10$. But the answer for this question set up a hypothesis of $H_0: \mu = 10,H_a: \mu>10$. Shouldn't it be the first way?


title: Difference between fixed effects models in R (plm) and Stata (xtreg)

text: I'm trying to re-create an analysis done using `Stata` function `xtreg` (though I don't have the code) with `R` package `plm`, and I'm having trouble translating between the two. A minimal example of the model I'm trying to estimate would be as follows: $$ \text{Dependent}_{i,t} = \beta_0 + \gamma_i + \delta_t + \beta_1*\text{Dependent}_{i,t-1} + \beta_2*\text{Independent}_{i,t} + \epsilon_{i,t} $$ where $\gamma_i$ is a control for the cross-sectional factor (States), $\delta_t$ is a control for each time period (Year), $\text{Dependent}_{i,t}$ is a dependent variable in the form of a proportion, and $\text{Independent}_{i,t}$ is an attribute of a record that changes with both the cross-sectional and time component of the data. I've attempted to estimate the coefficients for the model using `R` with the following (generalized) code:               model <- plm(Dependent_it ~ Independent_it + lag(Dependent_it,1),                   data   = data,                   method = "within", #fixed effects model                  effect = "twoway", #does the gamma_i and delta_t parts (I think)                  index  = c("State", "Year")                  )      I've tried to re-create it in `Stata` using the following (after running `xtset State Year, yearly` of course)               xtreg Dependent Independent l.Dependent i.Year, fe      I get different results, so I must be doing something wrong with either `xtreg` or `plm`, or both. Can someone enlighten me as to if I've coded the above model appropriately in either language?


title: Significance of combined lagged and contemporaneous variables

text: In several analyes I performed, I occasionally observed that when dealing with lagged variables, including contemporaneous values of the same variables in the regression increases the signficance of the lagged variables. The contemporaneous values are insignificant themselves nevertheless. For example, I regress market returns on the lagged value of a variable, i.e. X(t=-1). The coefficient is significant and also robust to other controllvariables. Nevertheless, when I include the value X(t=0), the coefficient gains substantial significance, e.g. p jumps from 2.4% to 0.1%. The coefficient on X(t=0) is insignificant though, but adj. R² increases. I am not sure what to make out of this. Should I include X(t=0) as long as there is economic intuition to it? Best Micha


title: Discrete features are less important

text: I am trying to train a random forest classifier. As predictors, I keep both discrete features and continuous features (the discrete ones including booleans, counters, etc., and the continuous contains floats). Training the model both in `R` and in Python's `scikit.learn`, I get that the discrete features are always less important (using the importance feature of the model). I wonder if that's a characteristic of my data or perhaps there's some innate characteristic of discrete features. I am not even sure it makes sense to mix discrete and continuous features in the same predictors matrix. Some numbers: a current test with about 30 features (not much, it'll grow larger) yields mean importance 0.035; the continuous features importances range from 0.037 to 0.06 (mean > 0.04) and the discrete features importances range from 0.005 to 0.03 (mean < 0.02). To sum, I have two questions:   1. is there something "inherent" in discrete features that makes them "less important" (according to the algorithm for importance calculation)   2. does it make sense to mix the feature types?


title: How do I specify a loss matrix in rpart?

text: I have a basic question - how to I set up a loss matrix to weigh the cost of a false positive higher than a false negative? I am trying to produce a tree in rpart to classify a disease with high specificity. Any advice would be appreciated, thanks


title: Confidence region about a regression model

text: I have a set of means and standard deviations. For each mean I can calculate a 95% confidence interval. I plot these means and confidence intervals against an independent variable and I fit a best sigmoidal curve through this data. I want to calculate the region of 95% confidence above and below this curve (i.e. be able to make estimates of the 95% confidence interval between the measured data points.) What should I be looking to use? Should I fit a second sigmodal curve to the upper limits of each data point's interval and a third curve to the lower limits of each interval?


title: Monte Carlo integration aim for maximum variance

text: I have a question about Monte Carlo integration. As I understand it the method takes a region S of known volume V which contains the region T specified in the definite integral. $T \in S$. Then random points in S are selected and it is checked if they belong to the volume that is the volume under the graph. Whether a point falls in the volume or not is a binomial random variable and if 30% of the sampled points do lie in the volume under the graph then the estimated volume is 0.3V. What I have read is that the optimize the method it is best to define S so that there is the greatest variance when sampling, i.e. when the probability of a point lying in the volume under the graph is 0.5. This seems counter-intuitive as the greatest variation would cause the greatest sampling error and the widest confidence interval. However when I tried to estimate the integral $\int_{x=2}^{3} \int_{y=6}^{7} xy\ dy\ dx$ I tried it in a region where the estimated volume was about 50% of the volume of S and in a region where the estimated volume was about 2% of S. I repeated the simulation several times and I was getting more accurate results with the 50% version. I am stumped as to why this happened and why it is recommended to use a figure of around 50%. Can someone please explain why having a large variation is advantageous.


title: Confused with MCMC Metropolis-Hastings variations: Random-Walk, Non-Random-Walk, Independent, Metropolis

text: Over the past few weeks I have been trying to understand MCMC and the Metropolis-Hastings algorithm(s). Every time I think I understand it I realise that I am wrong. Most of the code examples I find on-line implement something that is not consistent with the description. i.e.: They say they implement Metropolis-Hastings but they actually implement random-walk metropolis. Others (almost always) silently skip the implementation of the Hastings correction ratio because they are using a symmetric proposal distribution. Actually, I haven't found a single simple example that calculates the ratio so far. That makes me even more confused. Can someone give me code examples (in any language) of the following:   * Vanilla Non-Random Walk Metropolis-Hastings Algorithm with Hastings correction ratio calculation (even if this will end up being 1 when using a symmetric proposal distribution).   * Vanilla Random Walk Metropolis-Hastings algorithm.   * Vanilla Independent Metropolis-Hastings algorithm. No need to provide the Metropolis algorithms because if I am not mistaken the only difference between Metropolis and Metropolis-Hastings is that the first ones always sample from a symmetric distribution and thus they don't have the Hastings correction ratio. No need to give detailed explanation of the algorithms. I do understand the basics but I am kinda confused with all the different names for the different variations of the Metropolis-Hastings algorithm but also with how you practically implement the Hastings correction ratio on the Vanilla non-random-walk MH. Please don't copy paste links that partially answer my questions because most likely I have already seen them. Those links led me to this confusion. Thank you.


title: Interpreting the random effect in a mixed-effect model

text: I am looking at several dependant variables for which I created LMMs of the following kind:               DV ~ Group + (1|Subject) + (1|Time)      Now I am struggling with how to interpret the output concerning the random effects, e.g. for DV1:               Random effects:      Groups   Name        Variance Std.Dev.      Subject   (Intercept)  6.2     2.49          Time      (Intercept) 67.7     8.23          Residual             50.2     7.08          and for DV2:                Groups   Name        Variance Std.Dev.      Subject   (Intercept)  27.2     5.22         Time      (Intercept)  17.7     4.21         Residual             125.2    11.19         Are conclusions such as the following valid?   1. Variability between subjects was about twice as high for DV2 than DV1.   2. Variability between dates (time) was higher for DV1 than DV2.   3. For DV1 time accounted for a higher variability than subject   4. For DV2 a similar amount of variability can be attributed to time and subject.   5. For DV1 there was a variability of about 8 [units] between years.   6. For DV2 variability between subjects was about 5 [units]. Thanks!


title: Heteroskedasticity, autocorrelation robust standard errors for SPSS

text: Is there a way of performing HAC robust standard errors in SPSS?


title: SSCP in analysis of covariance

text: I want to use analysis of covariance in R. The model is $y_{ij} = \mu + t_i + r_j + b(x_{ij} - \bar{x}) + e_{ij} $ In R it is possible to get Ryy, Tyy, Eyy, SSE as well as b and adjusted treatment or block with using 'lm' and 'anova'. However I need Rxx, Txx, Exx, Rxy, Txy and Exy (SSCP). SAS returns these values. Is there any way to obtain these values in R. Thanks.


title: Fixing error variance in a path analysis to model measurement error in scales using sem package

text: I want to construct a path analysis model that can account for measurement error in totally aggregated parcels, which refer to parcels where all of the items in a scale are summed or averaged. If I am not mistaken, Bollen (1989) advocates the following formula for computing the error variance of each parcel: $$\big(1 - \alpha(\text{parcel})\big)\times \text{Var}(\text{parcel}),$$ such that $\alpha$ refers to Cronbach's alpha, which is a measure of reliability. What follows is the following path analysis model in theory:               path.inf.final <- specifyModel()     pRU -> sRU, test1     pRU -> rRU, test2     sRU -> rRU, test3     sRU -> power_alt, gam1     pRU -> power_alt, gam2     rRU -> power_alt, gam3     sRU -> ms_alt, gam7     pRU -> ms_alt, gam8     rRU -> ms_alt, gam9     sRU <-> sRU,  NA, (1 - alpha(sRU))*(variance(sRU))     pRU <-> pRU, NA, (1 - alpha(pRU))*(variance(pRU))     rRU <-> rRU, NA, (1 - alpha(rRU))*(variance(rRU))     power_alt -> em, iota1     power_alt <-> power_alt, NA, (1 - alpha(power_alt))*(variance(power_alt))     ms_alt -> em, zeta1     ms_alt <-> ms_alt, NA, (1 - alpha(ms_alt))*(variance(ms_alt))     em -> relM1, eta1     em <-> em, NA, (1 - alpha(em))*(variance(em))     relM1 <-> relM1, lam1      After testing this model, the model fits were unnaceptable and package sem could not compute modification instances. In contrast, when I do not fix error variances, after a few modifications to the model, I am able to reach a theoretical cogent model that has acceptable model fits. I have several questions on this regard:   1. Am I understanding Bollen (1989) correctly and using the formula correctly in my path analysis model?    2. I also heard that rather than Bollen's formula, I could use alpha or alpha squared. Is this correct? If so, would it follow that I just input it into the model in much the same manner like I used Bollen's formula. I also tried this approach and reached unsatisfactory results.   3. The dependent variable in the path analysis model is binary and I could argue that there is no measurement error associated with it. How could I fix the error variance for this variable, i.e., relM1, such that the aforementioned is accounted for. Any help would be much appreciated. References:   * Bollen, K. A. (1989). _Structural equations with latent variables_. New York: Wiley.


title: Bayesian plots: construction and interpretation

text: I've heard a couple of definitions of a Bayesian plot, but I am not 100% sure what it is and whether it is possible to plot this in an R. An example of the plot is here: http://www.springerimages.com/Images/MedicineAndPublicHealth/1-10.1007_s11307-008-0154-3-1 In this particular example there are four different tests, the prior probability is evaluated based on the proven probability of cancer recurrence based on 7 patients. This is compared against varying prevalences of tumour (As far as I can understand this is the the x-axis reflects these changing prevalences). The y-axis is the posterior probability. What I don't fully understand is what the likelihood is and how they get it from the different methods to actually calcualte the posterior and show there are differences between these methods. I am not necessarily interested in understanding this particular example. I just want to understand how I can construct and interpret this plot on my own data. I understand that my prior should be my disease prevalence, but what I don't understand is how the likelihood should be represented in this plot (is it the x-axis?). From my own data I can evaluate the posterior because I know the prior and the likelihood... I just don't know how to put them into this plot. Has anyone made such a plot for their own work that can help me to understand the intuition?


title: Adaptive Regularization of Weights (AROW): classification dependent on feature index

text: I am trying to classify texts using a bag of words algorithm. My feature vector is thus a large array of words. In order to build my feature vector, I parse the sample and put a one in the index corresponding to the word found. I use sparse format to represent the feature vectors. A little example to clarify:               feature vector = ('cat', 'dog', 'parrot', 'cow', 'fish')     sample = 'The cat is playing with the cow'     vector =  [0:1, 3:1]      I am using two classifiers in my tests: AROW and SVM. I found something very strange with the arow classifier: if I change the index of the words in my feature vector, classification results vary enormously. For example, if I change my initial feature vector to be:               feature vector = ('fish', 'parrot', 'cat', 'dog', 'cow')      My results won't be the same when I am using AROW. On the other hand, SVM gives exactly the same results. I am using the freely available AROW library http://code.google.com/p/arowpp/ in my tests. Is this a "feature" of AROW or is it likely an implementation bug in the open source library?


title: What is the probability of rolling all faces of a die after n number of rolls

text: It is fairly easy to figure out what is the average number of rolls it would take to roll all faces of a die [$1 + 6/4 + 6/4 + 6/3 + 6/2 + 6/1 = 14.7$], but that got me thinking of a seemingly more complicated problem. If you roll a die 1-5 times, the odds of _all_ faces showing is obviously 0. If you roll a die 6 times, the odds of all faces showing can easily be calculated as $$1 * (5/6) * (4/6) * (3/6) * (2/6) * (1/6) = .0154$$ or 1.54% This is where I get stuck. How to do it $n=7$ or more times and calculate it for general $n$? Any tips are helpful!


title: Value for which the PDF(value) is maximal in a distribution with skew?

text: I am working on a project where I need to chart statistical data and related, skewed distributions a la http://en.wikipedia.org/wiki/Skew_normal_distribution. Unlike with normal distributions, in these charts, when there is skew, neither the mean nor the median appear to be at the maximal point of the probability distribution curve drawn along the graph's domain, which is fine, but does this point have a name? And more importantly, is there a more direct way to discover the value along the domain of values such that the PDF(value) has the maximal value? I can currently find this point with a ~~binary search~~ hill-climbing algorithim, but, this seems a rather round-about way to go about it. While this point may or may not have any intrinsic value in statistics, from the standpoint of scaling the resulting graphs, I'd like to calculate the highest point of the distribution. I am not a statistician, so it could very well be that the answer is on the very Wikipedia page I posted above, but I am not familiar with the parlance of this level of Statistics. **Summary** The consensus is that there is no closed form of computing the mode from an arbitrary Skew Normal Distribution, so, a Hill Climbing algorithm is probably a sensible approach. A good starting point for the algorithm, seems to be the distribution's median, which does have a closed form, and except in cases of extreme skew, is fairly close to the mode. At extreme skews, the location is probably the best starting point. I have created a nifty animation in R to visualise this: ![enter image description here](http://i.stack.imgur.com/VBvys.gif) PS: R is fun! =)


title: How to use one variable in regression with many independent variables in lm()

text: I need to reproduce this code using all of these variables.               composite <- read.csv("file.csv", header = T, stringsAsFactors = FALSE)     composite <- subset(composite, select = -Date)     model1 <- lm(indepvariable ~., data = composite, na.action = na.exclude)     composite is a data frame with 82 variables.      UPDATE: What I have done is found a way to create an object that contains only the significantly correlated variables, to narrow the number of independent variables down. I have a variable now: sigvars, which is the names of an object that sorted a correlation matrix and picked out only the variables with correlation coefficients >0.5 and <-0.5. Here is the code:               sortedcor <- sort(cor(composite)[,1])     regvar = NULL          k = 1     for(i in 1:length(sortedcor)){       if(sortedcor[i] > .5 | sortedcor[i] < -.5){         regvar[k] = i       k = k+1      }     }     regvar          sigvars <- names(sortedcor[regvar])      However, it is not working in my lm() function:               model1 <- lm(data.matrix(composite[1]) ~ sigvars, data = composite)      Error: Error in model.frame.default(formula = data.matrix(composite[1]) ~ sigvars, : variable lengths differ (found for 'sigvars')


title: How to use SVD for dimensionality reduction to reduce columns specifically?

text: My original data has many more columns (features) than rows (users). I'm trying to reduce the features of my SVD (I need all of the rows). I found one method of doing so in a book called "Machine Learning in Action" but I don't think it will work for the data I am using. The method is as follows. Define SVD as A = USV^T. Set an optimization threshold (i.e., 90%). Calculate the total sum of the squares of the diagonal sigma matrix. Calculate how many sigma values it takes to reach 90% of the total sum of squares. So if that turns out to be 100 sigma values, then I would take the first 100 columns of the 'U' matrix, first 100 rows of the 'V^T' matrix, and a 100 x 100 square matrix out of the 'S' (or Sigma) matrix. I would then calculate A = USV^T using the reduced matrices. However, this method targets the rows of my original data, since the dimensions of the Sigma matrix are rows x rows. How would I target the columns of my original matrix?


title: Wilcoxon test to support logistic regression

text: I want to test the hypothesis that IV 2 (post-test value) is a better predictor for the criteria (no success vs. success) than IV 1 (pre-test value). Those IVs are indicators for static vs. dynamic learning. I conducted a logistic regression analysis for both IVs seperately and it seems that IV 2 is indeed the better predictor. Now I would like to support this finding with the wilcoxon-test. I seperated the data into the groups "no success" and "success", but I'm not quite sure how to go from here because I have trouble to understand which procedure is correct to test my hypothesis... Do I have to compute a variable with a hypothetic value of 0 ---> wilcoxon test with IV 1 and the hypothetic value variable, repeat this step with the IV 2 Or do I have to follow the usual procedure and just compare IV 1 and IV 2? What is the right procedure to complement the logistic regression analysis? Thanks for any advice!


title: Likert Scale Analysis - Pre/Post

text: I was hoping somebody could help with determining the correct statistical test to use. Basically a teaching session was done which assessed confidence and perceptions of handling a situation using 3 questions. Each question had a Likert Scale (1-5). The same 3 questions on the pre and post questionnaire. The problem is that not everybody filled out the pre/post questionnaire so the numbers are different and potentially more in one group than the other. Is there any way to analyse this sort of data


title: How to test my data against an specific normal distribution?

text: I need to test my data to see if it follows a normal distribution with specific mean and std like N~(mu, std) I know that this can be done by Kolmogorov-Smirnov test which has a function in both MATLAB and R but the default for these function is standard normal and I do not know how should I specify my mu and std in the function. If you can help me in any of these two softwares I will be thankful.


title: Independent paired t-test -> Mann Whitney U

text: Just a quick introduction, I'm rubbish at stats. I asked a question before which was very well answered but I have since modified the experiment in a way that confuses me as to which test I should use. So in the previous question we finished concluding that an independent t-test would be suffice as I am comparing the effect of irradiation on each leaf so, one half would be the control or the "before", and the irradiated half would be the "after" treatment. Well, I got my results, but I didn't map out which irradiated half went with what control as it was very difficult to keep track so I ended up with a jumbled set of controls to a jumbled set of irradiated. So I made the following changes to my experiment: 1) Instead of having 5 circles per half of leaf, I instead opted for just 1. This meant it was quicker to get results as things like overlapping circles reducing light intake for the leaf discs on bottom didn't happen. So I would initially cut all my leaves in half. take one sample from one half, and whilst waiting for photosynthesis to take place I would irradiate the other halve(s), which brings me on to my second change.. 2) I did 3 syringes at the same time instead of just 1 quadrupling my speed. (That's what I meant when I said I did four at the same time), so 4 control, whilst waiting, irradiate other halves of leaves for a minute and switch once photosynthesis finished, then switch again to another batch of 4 leaves and repeat. I ended up with 16 results either side


title: Interpretation of fixed effect model (panel data)

text: I am doing a work to school in which I have to analyze effect of expenditure for research and development on Gini index. I have _Panel data_ and when I do _Fixed effect model_ , I have problem with interpreting them. If anybody can help, I am looking for some place where I can find what each of the results I get means because I have to do other models as well. This is what I get from fixed effect model:               Oneway (individual) effect Within Model          Call:     plm(formula = log(Gini) ~ log(VaR), data = data, model = "within")          Balanced Panel: n=17, T=11, N=187          Residuals :        Min. 1st Qu.  Median 3rd Qu.    Max.      -0.3140 -0.0933  0.0149  0.0972  0.2270           Coefficients :                Estimate Std. Error t-value  Pr(>|t|)         log(VaR) -0.0272357  0.0078724 -3.4596 0.0006847 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Total Sum of Squares:    2.8467     Residual Sum of Squares: 2.6584     R-Squared      :  0.066139            Adj. R-Squared :  0.059773      F-statistic: 11.9691 on 1 and 169 DF, p-value: 0.00068471


title: Linear model for this problem? Which regression coefficients should I use?

text: > Suppose a researcher was investigating gender-based wage disparities across > 10 different firms. The researcher obtained the following data from these 10 > companies employees: 1) years of experience; 2) education; 3) gender; and 4) > the employees’ company. Suppose the researcher is interested in whether > there are overall gender-based wage differences across the 10 companies > after controlling for experience and education. Describe the model you would > construct to test this hypothesis and provide an interpretation of > parameters of interest. I was thinking that I would use something like this: Fit1: $Wage = B_{0} + B_{1}Gender + B_{2}YearsofExperience + B_{3}Company$ Fit2: $Wage = B_{0} + B_{2}YearsofExperience + B_{3}Company$ Test using F = Anova(Fit1,Fit2). $H_{0}$: Wage is the same, reject if $F > F^*_{.05}(1,4)$ One question that I have about the above: Is doing $B_{3}Company$ correct? Or do I need to have a B value for each company, ie $B_{3}Company1, B_{4}Company2, B_{5}Company3$ ... etc. Also: > Suppose a researcher was investigating gender-based wage disparities across > 10 different firms. The researcher obtained the following data from these 10 > companies employees: 1) years of experience; 2) education; 3) gender; and 4) > the employees’ company. Suppose the researcher is interested in whether > there is evidence that the wage gap differs across companies after > controlling for experience and education. Describe the model you would > construct to test this hypothesis and provide an interpretation of > parameters of interest. For this problem, how would I go about making the predicted value "wage difference" versus just "wage" like before?


title: Test of whole distribution for two samples

text: I have a fundamental question: is it reasonable to use a statistical test to find out if 2 datasets are similar or not? Some comments have called out the T-test as being unable to answer that question. Why is the T-test insufficient? Furthermore, what other tests are there to determine differences in distributions for the two sample problem.


title: statistical test for sequence of data

text: There is a sequence of data. It might have different types of trend, i.e., increasing, decreasing, and even cycling pattern. But there is an assumption that these trends should be smooth. (maybe not mathematically strict). In other words, there should not have a spike, or a sharp changing point; or the changing ratio $\frac{Y_{N+2}-Y_{N+1}}{Y_{N+1}-Y_N}$ should not change sharply. Are there any statistical test that can define( or detect) this kind of smooth pattern.


title: Can all methods made for certain level in measurement hierarchy used at all higher levels?

text: Let's say there's a measurement hierarchy for discrete variables, that looks like nominal (grade pass/fail) < ordinal (grade A/B/C) < interval (grade 90-100; 80-89; 70-79, etc). Hierarchy is used to refer to the possiblity that we can assign order to the discrete values. I read a line that methods applicable for one type of variable can be used for the variables at higher levels too, but not at lower levels. Can any one explain and exemplify why this is true (if it is)?


title: How to create a weekly alert system that works for different sample sizes?

text: I'm trying to build a weekly alert system to let me know when something unusual has happened in my analytics data that might require further investigation. Currently there is an outdated system in place that only looks at week over week % increase, which obviously is not working. Let's think about the only metric as website hits (x). I have about 2,000 websites that show great difference in size from each other. We can categorize them under 2 behaviors.   1. Shows high activity with low variance (x~ 100,000)   2. Shows low activity with high variance (x~ 100) I am trying to come up with an approach that is robust for all and minimize false positives. Need something smart enough to tell me that 10% increase in 100K is more important than 50% increase in 100.   * What do we have in literature for these kind of problems? What are they called?    * Is it better to categorize the websites based on their activity (average, std_dev) and create separate alert system?   * What would the best way be to deal with seasonality? (ie travel seasons which would skew averages and standard deviations)   * How should I deal with small websites? Would I have to take a longer time frame to calculate averages? How would I decide how many to take?


title: Predicting a Certain Type of failure and deciding the input time series

text: I am trying to predict the time to certain type of failure given the following data on Certain Factory Equipments. The data I have are readings collected every day for sensor installed on those equipments .On Same day, an equipment can have different Repaires performed,and on some days no repairs. When a repair is performed it means some part of that equipment failed. Our goal is to find the time to failure of that part so we can perform preventative maintenance.I converted the data so I have only one type of failure in my data , Part1 failed or not failed is my binary output. I intend to do sperate analysis for diffferent part failures. I was thinking of using the randomForestSRC package to fit a model.I have not done any survival Timeseries Modeling. my challenge is how to define the Input Time series in the surv(time,output) What other things should be kept in mind when doing survival analysis on such data?? Can i just sort the data by time and use that as input or I need to define my Starttime (time when the equipment was brandnew and stop time as time when the equipment had certain part failure??) Should the analysis be done separately on each Equipment ID ?? Any ideas appreciated. Sample data . ![enter image description here](http://i.stack.imgur.com/dfTNM.png)


title: Conservative confidence interval

text: One month before the election, a poll of a large number of randomly selected voters showed 65% planning to vote for a certain candidate. The newspaper article reported a 90% conservative margin of error of 4 percentage points. Use the above information to construct a 90% conservative confidence interval for the population proportion of all voters planning to vote for this candidate. Can someone please explain why the answer is 0.65 +/- 0.04?


title: Should the predictor variables be normally distributed for Poisson glm?

text: This is probably a really basic question, but it's the first time I've created a model that defines Poisson as its error family. In setting up my variables to make the model, should I be concerned about whether or not predictor variables are normally distributed, and if not, should I be attempting to transform them to make them normal? Or alternatively, should the residuals of simple regressions between each predictor and the response variable be normally distributed? Or is this something I look at overall, once the model is made, by looking at a histogram of the residuals of the full model? Or, does the normal distribution not even apply in this case because I have specified Poisson errors? Thanks!


title: Is this a case for ordered logistic regression?

text: I have the following study setup: Three groups of people were asked a question, and the answer was ordinal (likely, somewhat likely, somewhat unlikely, unlikely). In my data set, I have a contingency table between the response and groups, and would like to know whether some of the groups of people were more inclined to answer the question with likely than others. EDIT: I didn't clearly state this above, but I want to treat the answer as ordinal, i.e. likely = 1 < somewhat likely < somewhat unlikely < unlikely. I think I could do it with an ordered logistic regression model, but isn't that like shooting cannon at sparrows?


title: Test statistic as Gamma distribution

text: I want to perform a simulation study where I want to use gamma random variable as the test statistic for some test. Does anyone know a hypothesis test whose test statistic has a gamma distribution. Thank you very much.


title: Histogram for a sample in Python easily?

text: Instead of this simplified-version of this file:               #!/usr/bin/env python     import numpy as np     import pylab as P               data = [x.rstrip().split("\t") for x in file(".data", "r").readlines()]     sdData = [float(x[1]) for x in data]     retData =[float(x[0]) for x in data]          # Calculating medians for estimating mean and sd     sdMed = sorted(sdData)[len(sdData)/2]     retMed= sorted(retData)[len(sdData)/2]               mu, sigma = retMed, sdMed     x = mu + sigma*P.randn(10000)          # the histogram of the data with histtype='step'     n, bins, patches = P.hist(x, 20, normed=1, histtype='stepfilled')     P.setp(patches, 'facecolor', 'g', 'alpha', 0.75)          # add a line showing the expected distribution     y = P.normpdf( bins, mu, sigma)     l = P.plot(bins, y, 'k--', linewidth=1.5)               #     # create a histogram by providing the bin edges (unequally spaced)     #     P.figure()     P.show()      is there something like this `hist(SampleDataPoints)`?


title: Logistic predict.glm() function error: factor has new levels

text: I fit a logistic on three numeric continuous variables, followed by a categorical factor [Y, N].               logit2A <- glm(DisclosedDriver ~ VehDrvr_Dif+POL_SEQ_NUM+PRMTOTAL+SAFE_DRVR_PLEDGE_FLG, data = DF, family = "binomial")       Fit looks wonderful.               Coefficients:                         Estimate Std. Error z value Pr(>|z|)         (Intercept)           -2.204e+00  2.253e-01  -9.782  < 2e-16 ***     VehDrvr_Dif            2.918e-01  1.026e-01   2.845 0.004440 **      POL_SEQ_NUM           -1.893e-01  5.617e-02  -3.370 0.000751 ***     PRMTOTAL               1.109e-04  5.526e-05   2.006 0.044804 *       SAFE_DRVR_PLEDGE_FLGY -7.220e-01  1.633e-01  -4.422 9.76e-06 ***      So obviously R took the Safe_Drvr_Pledge_Flg categorical factor variable and placed all 'N' values in reference or intercept as opposed to the listed 'Y'. Now I want to take my fit and calculate the probabilities that my model determines. And here comes the error:               > DF$P_GLM<- predict.glm(logit2A, DF, type="response", se.fit=FALSE)         Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels) :      factor SAFE_DRVR_PLEDGE_FLG has new levels       Umm... no it doesn't, because I just fit the model with the exact same data I'm trying to use for the prediction. What's the problem? Trying to respond to first comment: Don't know what you mean. I've got 3500 rows of data... It's a logistic regression on 4 continuous variables and one categorical. The categorical has two values, Y or N. My glm fit give the numbers given. I just want to plug it all back in with the predict function and it gives me that error. Here's the categorical variable:                > DF$SAFE_DRVR_PLEDGE_FLG      [1] Y Y Y Y Y Y Y N Y Y Y Y Y Y Y Y Y Y Y N Y Y N Y Y N Y Y N Y Y Y Y Y Y Y Y Y Y N Y Y N Y Y Y N Y Y Y Y Y N Y Y Y Y Y Y      [60] Y Y Y Y N Y Y Y Y Y Y Y Y N Y Y Y N N Y N Y Y Y Y Y N Y Y N Y N N Y Y Y N Y Y Y Y N Y Y Y Y Y N Y N Y N Y Y Y Y Y N Y      [119] N Y Y Y Y Y Y Y Y N Y Y Y Y Y Y N Y Y Y N Y Y Y N Y Y Y N N Y N N N Y N Y Y Y N N Y Y N Y Y Y Y N N Y Y Y Y N N Y N N      Levels:  N Y      What do you mean by a working example? The fit works. The probability output of the predict function doesn't...


title: Confusion related to Gibbs sampling

text: I came across this article where it says that in Gibbs sampling every sample is accepted. I am a bit confused. How come if every sample it accepted it converges to a stationary distribution. In general Metropolis Algorithm we accept as min(1, p(x*)/p(x)) where x* is the sample point. I assume that x* points us to a position where the density is high so we are moving to the target distribution. Hence I suppose that it moves to the target distribution after a burn in period. However, in Gibbs sampling we accept everything so even though it may take us to a different place, how can we say that it converges to the stationary/target distribution Suppose we have a distribution $p(\theta) = c(\theta)/Z$. We cannot calculate Z. In metropolis algorithm we use the term $c(\theta^{new})/c(\theta^{old})$ to incorporate the distribution $c(\theta)$ plus the normalizing constant Z cancels out. So it's fine But in Gibbs sampling where are we using the distribution $c(\theta)$ For eg in the paper http://books.nips.cc/papers/files/nips25/NIPS2012_0921.pdf its given so we don't have the exact conditional distribution to sample from, we just have something that is directly proportional to the conditional distribution ![enter image description here](http://i.stack.imgur.com/8yH9o.png)


title: How to visualize what canonical correlation analysis does (in comparison to what principal component analysis does)?

text: Canonical correlation analysis (CCA) is a technique related to principal component analysis (PCA). While it is easy to teach PCA or linear regression using a scatter plot (see a few thousand examples on google image search), I have not seen a similar intuitive two-dimensional example for CCA. How to explain visually what linear CCA does? Your Suggestions? EDIT: rephrased the question a bit (see comments)


title: Significance testing for different groups in multiple regression with dummy-coded interaction variables

text: Consider a multiple regression predicting outcome `Y` using a continuous predictor `P` and an interaction between `P` and a group variable `G` with the values 0, 1, and 2. We create two interaction variables, `PxG1` (set to P when G is 1 and 0 otherwise) and `PxG2` (set to P when G is 2 and 0 otherwise). When P, PxG1, and PxG2 are used as simultaneous predictors of `Y` (assume we're not interested in the effect of G by itself), the regression coefficient for `P` will give the relationship between P and Y for group 0, and the coefficients for `PxG1 and`PxG2` will give the difference between the relationship for group 0 and the relationship for groups 1 and 2, respectively. Question 1: The significance test for the coefficient of P can be interpreted as showing the significance of the correlation between P and Y in group 0, right? Question 2: What if I want to find out whether the correlation between P and Y is significantly different from 0 in the other two groups? Could I switch to using the interaction variables `PxG0` and `PxG2`, making G1 the reference group, and interpret the coefficient of P as showing the significance of the correlation in group 1? Question 3: If I do that, does it create multiple comparison issues? I think the answer is no, because the information contained in P and any two of the interaction variables fully constrains the information in the third, so all I'm doing is parsing the available information in different ways. But it seems odd that rearranging my variables can provide additional significance tests.


title: Fast approximation to inverse Beta CDF

text: I am looking for a fast approximation to the inverse CDF of the Beta distribution. The approximation need not be precise, but more stress is on simplicity (I'm thinking Taylor expansion of the first 1 or 2 terms). Any ideas? A few more details: in my Beta distribution, $a,b$ are integers. Also the value I am evaluating the inverse of is always 0.05.


title: Higher order generalization of the multivariate normal distribution

text: In some sense the multivariate normal is the "nicest" distribution that we can describe using only a vector (rank one tensor) and a symmetric positive definite matrix (rank two tensor). $$\mathcal{N}(\mu_i, \Sigma_{ij}) = |2\pi\Sigma|^{-\frac{1}{2}}e^{-(x-\mu)^T\Sigma^{-1}(x-\mu)}$$ Is there a generalization to higher order tensors? I.e., is there a "nicest" distribution described by a vector, symmetric matrix, and symmetric third order tensor? $\mathcal{M}(\mu_i, \Sigma_{ij}, \Pi_{ijk})$ ?


title: How to succinctly report the correlation between two variables for 20 different subsets of the data in a table?

text: I have 20 correlation coefficients (r) to report, for example, .34, .45, .67, .23. and so on. These are correlations between Variable A and Variable B, based on gender, age, education levels etc. I am after a succient method to show all the correlations in one glance. Is there any best practice in regard to reporting all these correlations in a **table format**?


title: Poisson binomial distribution PMF conjecture

text: I have an intuition that re-casting the Poisson binomial distribution probability mass function in the following way might lead to a more elegant closed form than the DFT approach: Sum the $\log_2$'s of the probabilities of each independent trial. The negation of that sum is equivalent to a new number of trials each of which has a $p=1/2$. Of course, if all the $p$'s start as $1/2$ the number of trials will be unchanged and one can degenerate the existing formulas accordingly to yield an $O(n)$ solution. But if the $p$'s differ, something else must be done (for example, because not every $p=1/2$ trial is independent). My conjecture (or perhaps "wishful thinking" is more appropriate) is that there is an $O(n)$ closed form which will tell how one must re-number the desired number of successful trials ($k$) so it can be applied to the new number of trials where each has a $p$ of $1/2$, but not all are independent, yielding the same result as one would get from the ordinary PBD pmf.


title: R and SPSS differences in pca loadings

text: I performed two principal components analyses: in R and in SPSS - using the same dataset and the same variables. I got the same results - at least to some point. The eigenvalues are the same (I used the correlation matrix in both cases, no rotation), but the loadings are different. I decided to take a look at some plots, so below you can see two plots of those R and SPSS pca loadings. The first plot presents a little bit changed loadings of pca in R. Namely, these are the opposite signs of R loadings .   The second plot presents the original loading of pca in SPSS. These plots are the same, but I can't figure out why they different in numbers. Does R or SPSS do any implicit loadings transformation? ![enter image description here](http://i.stack.imgur.com/ZOugY.jpg)


title: $\alpha,\beta\ge 1$ in a Beta distribution. What does it imply for the mean and variance?

text: The Beta distribution has the PDF: $$f\left(x\right)=\frac{x^{\alpha-1}\left(1-x\right)^{\beta-1}}{\mathrm{B}\left(\alpha,\beta\right)}$$ for $0<x<1$, and $f(x)=0$ otherwise. The parameters $\alpha,\beta$ are positive real numbers. The mean and variance are given by: $$\mu=\frac{\alpha}{\alpha+\beta},\quad\sigma^{2}=\frac{\alpha\beta}{\left(\alpha+\beta\right)^{2}\left(1+\alpha+\beta\right)}$$ which can be inverted to give $\alpha,\beta$ in terms of the mean and the variance as $\alpha=\lambda\mu$ and $\beta=\lambda\left(1-\mu\right)$, where $$\lambda=\frac{\mu\left(1-\mu\right)}{\sigma^{2}}-1$$ Now I want to impose the condition that $\alpha,\beta \ge 1$. What does this imply for the mean and the variance? That is, is there a simple condition on $\mu,\sigma^2$ that is equivalent to $\alpha,\beta \ge 1$?


title: Age period cohort shift parameter

text: Consider the companion papers by Clayton & Schifflers (1987) here: http://www.pauldickman.com/cancerepi/handouts/handouts_time_trends/clayton_schifflers_1987a.pdf, http://www.pauldickman.com/cancerepi/handouts/handouts_time_trends/clayton_schifflers_1987b.pdf. Pages 464 and 474 discuss the drift parameter, but it's not clear how to compute the difference with the reference value using grouped data. Suppose my reference is 1960-64 and my period or cohort of interest is 1970-74. It seems like the difference in the midpoint, 10, is the proper choice here since no distinction is made between birth years falling within the intervals. Is that the value entered into the analysis for persons falling into these groups? Now, if one has individual level data, the exact differences could be computed, but would still need to be grouped in a manner similar to the rest of the data before conducting the analysis, correct?


title: Calculating statistically rare events from a set of observations

text: I have a Pandas data frame as follows: Event Location| Number of events | %age success A | 10 | 0.5 B | 2 | 1 C | 1 | 0 D | 100 | 0.3 E | 1 | 0 F | 1 | 1 ... ... The data represent the entire population of events. I am looking for a way to find Event Locations whose behavior is significantly different from others. Perhaps Detect samples with a high confidence of having a different mean than other Samples. To do this I am trying to calculate the probability of getting that sample in this distribution. For example, getting 1 event with either 0 or 1 success rate is likely. However, getting lots of events away from the mean is unlikely. I am thinking of using the t-test to find the statistical significant of that event sample in this population. Is there a better or another way to do this? To detect locations whose events are rare?


title: the effect of unbalanced classes

text: I have a data set with label "Y" and "N", which is splitted into training and testing sets. The ratio the labels in training and testing sets are 3:1 and 1:1, respectively. If I build model using the training data with random forest and 10-CV for model tuning, my question is, when the model is used predict the testing set, on average will the ratio of predicted labels be close to 1:1 or 3:1? I just want to understand how the unbalance in the training data influence the new data for prediction.


title: latent variables versus model parameters

text: I am quite confused with the distinction between a latent variable and model parameters. So say I have two observed variables $x$ and $y$ and they have some unknown relationship between them i.e. $y = f(x)$. Now based on that I have a generative model with $ y = f(x) + e$ where I model $e$ as zero mean independent and identically distributed noise with some variance $\sigma$. Am I correct in saying that: 1: The relationship $f$ that we try to estimate is the latent variable in this case. 2: $\sigma$ is a model parameter. Now, our estimate of $\sigma$ will also have some uncertainty associated with it. With these model parameters, are we saying that it has a fixed but unknown value that we try and estimate, but due to limited data size our estimates have uncertainty. However, underlying the model is the assumption that the value of a $\sigma$ takes a fixed value or is $\sigma$ also a random variable type entity?


title: How can I measure the association between 40 binary variables?

text: I have 40 binary (0/1) variables with 70 observations and would like to measure the association of occurrences between them. Is the presence or absence of a variable is associated with the presence / absence of another one? I believe I would be trying to do a multivariate analysis but with binary variables instead of continuous.


title: Interpretation of spatial lag models

text: I have estimated a linear regression model via R's `spdep` package (`lagsarlm` function), in which spatial auto-correlation was controlled for by an IDW matrix. I wonder whether the interpretation of model results is different from that of OLS models. For example, if a coefficient is 0.2, can I interpret it as: one unit increase in IV leads to 0.2 unit increase in DV?


title: How to calculate gelman and rubin convergence diagnostics in Excel

text: Does anyone know how to calculate the gelman and rubin convergence diagnostics in Excel. I'm trying to understand the steps in the formula, but not sure that I doing it correctly. I know it can be done in R, but they don't show the calculations. Thanks in advance. Ahmed


title: Some doubts about using GLMM

text: I'm analyzing the results of a hormone manipulation experiment. I measured a number of variables at three times in three groups. The groups are different sizes and not all individuals were measured every time, so I'm using GLMM rather than a repeated-measures ANOVA. I created the model then tested the significance of the terms (time, treatment, and time x treatment) with ANOVA. I'm quite new to GLMM, but after doing the tests, further reading suggests that my approach may be inappropriate, particularly with small data sets (I have ~seven animals per group). It seems that there is disagreement about what the degrees of freedom should be. This leads me to three questions: 1) Is this method acceptable? 2) If so, what would be an appropriate method to do post-hoc analyses to determine which groups differ? 3) Like the fake data, I have a number of negative results. Specifically, I often see significant time effects, but no effect of treatment or treatment X time. If I stick with this method, how can I calculate effects sizes and/or confidence intervals for such tests? Here are some fake data:               library(nlme)          datums<-data.frame(id=rep(1:20,each=3),var1=runif(60,4,6),var2=runif(60,25,30),var3=runif(60,0,1),var4=runif(60,10,15),var.time=rep(1:3,times=20),var.treatment=rep(c('a','b','c'),each=20))          datums$var.time<-as.factor(datums$var.time)     datums$id<-as.factor(datums$id)               #and now the GLMMs on each variable - I'll show just two here          var1.glmm<-lme(var1~var.time + var.treatment + var.time*var.treatment, data=datums, random = ~1| id)          var2.glmm<-lme(var2~var.time + var.treatment + var.time*var.treatment, data=datums, random = ~1|id)          summary(var1.glmm)          anova(var1.glmm)      I'm aware that the place for me to go is probably the Pinheiro and Bates book, but I don't have access to it at this time. Thanks in advance for any advice.


title: Multiple regression with no origin and mix of directly entered and stepwise entered variables using R

text: I am running a regression equation and I want to enter in 12 indepdendent variables then stepwise enter 7 more independent variables and not have an origin.   * DV is `shfl`.    * I want to enter in the following 12 independent dummy variables `ajan bfeb cmar dapr emay fjun gjul haug isep joct knov ldec`   * And then I want to enter in a stepwise fashion `slag6 slag7 slag8 slag9 slag10 slag11 slag12`   * And finally, I want there to be no origin. I've done simple regression, but nothing quite like this that enters in the primary variables and step enters several more.   * How can such a model be specified using R?


title: Why do the R functions 'PCA' and 'prcomp' give different eigenvalues?

text: You can use the decathlon dataset {FactoMineR} to reproduce this. The question is why the computed eigenvalues differ from those of the covariance matrix               > library(FactoMineR);data(decathlon)     > pr <- princomp(decathlon[1:10], cor=F)     > pr$sd^2           Comp.1       Comp.2       Comp.3       Comp.4       Comp.5       Comp.6      1.348073e+02 2.293556e+01 9.747263e+00 1.117215e+00 3.477705e-01 1.326819e-01            Comp.7       Comp.8       Comp.9      Comp.10      6.208630e-02 4.938498e-02 2.504308e-02 4.908785e-03       Using another method               > res<-PCA(decathlon[1:10], scale.unit=FALSE, ncp=5, graph = FALSE)     > res$eig               eigenvalue percentage of variance cumulative percentage of variance     comp 1  1.348073e+02           79.659589641                          79.65959     comp 2  2.293556e+01           13.552956464                          93.21255     comp 3  9.747263e+00            5.759799777                          98.97235     comp 4  1.117215e+00            0.660178830                          99.63252     comp 5  3.477705e-01            0.205502637                          99.83803     comp 6  1.326819e-01            0.078403653                          99.91643     comp 7  6.208630e-02            0.036687700                          99.95312     comp 8  4.938498e-02            0.029182305                          99.98230     comp 9  2.504308e-02            0.014798320                          99.99710     comp 10 4.908785e-03            0.002900673                         100.00000      Can you explain to me why the previously computed eigenvalues differ from those ? (the eigenvectors are the same)               > eigen(cov(decathlon[1:10]))$values      [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01      [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03      Thank you! **EDIT** Using the alternative `prcomp` method..               > prc <- prcomp(decathlon[1:10])     > prc$sd^2      [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01      [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03      Why do PCA and princomp give different eigenvalues?


title: Is it possible to analytically integrate $x$ multiplied by the lognormal probability density function?

text: Firstly, by analytically integrate, I mean, is there an integration rule to solve this as opposed to numerical analyses (such as trapezoidal, Gauss- Legendre or Simpson's rules)? I have a function $\newcommand{\rd}{\mathrm{d}}f(x) = x g(x; \mu, \sigma)$ where $$ g(x; \mu, \sigma) = \frac{1}{\sigma x \sqrt{2\pi}} e^{-\frac{1}{2\sigma^2}(\log(x) - \mu)^2} $$ is the probability density function of a lognormal distribution with parameters $\mu$ and $\sigma$. Below, I'll abbreviate the notation to $g(x)$ and use $G(x)$ for the cumulative distribution function. I need to calculate the integral $$ \int_{a}^{b} f(x) \,\rd x \>. $$ Currently, I'm doing this with numerical integration using the Gauss-Legendre method. Because I need to run this a large number of times, performance is important. Before I look into optimizing the numerical analyses/other pieces, I would like to know if there are any integration rules to solve this. I tried applying the integration-by-parts rule, and I got to this, where I'm stuck again,   1. $\int u \,\mathrm{d}v = u v - \int v \mathrm{d}u$.   2. $u=x \implies \rd u = \rd x$   3. $\rd v = g(x) \rd x \implies v = G(x)$   4. $u v - \int v \rd x = x G(x) - \int G(x) \rd x$ I'm stuck, as I can't evaluate the $\int G(x) \rd x$. This is for a software package I'm building.


title: STL trend of time series using R

text:
I am new to R and to time series analysis. I am trying to find the trend of a long (40 years) daily temperature time series and tried to different approximations. First one is just a simple linear regression and second one is Seasonal Decomposition of Time Series by Loess. In the latter it appears that the seasonal component is greater than the trend. But, how do I quantify the trend? I would like just a number telling how strong is that trend.                    Call:  stl(x = tsdata, s.window = "periodic")          Time.series components:             seasonal                trend            remainder                    Min.   :-8.482470191   Min.   :20.76670   Min.   :-11.863290365           1st Qu.:-5.799037090   1st Qu.:22.17939   1st Qu.: -1.661246674      Median :-0.756729578   Median :22.56694   Median :  0.026579468           Mean   :-0.005442784   Mean   :22.53063   Mean   : -0.003716813      3rd Qu.:5.695720249    3rd Qu.:22.91756   3rd Qu.:  1.700826647         Max.   :9.919315613    Max.   :24.98834   Max.   : 12.305103891              IQR:              STL.seasonal STL.trend STL.remainder data                 11.4948       0.7382    3.3621       10.8051            % 106.4          6.8      31.1         100.0            Weights: all == 1          Other components: List of 5        $ win  : Named num [1:3] 153411 549 365
         $ deg  : Named int [1:3] 0 1 1        $ jump : Named num [1:3] 15342 55 37
         $ inner: int 2       $ outer: int 0      ![enter image description here](http://i.stack.imgur.com/jwCSr.png)


title: How to initialize EM-algorithm when trying to fit data to a separable mixture model?

text: As a part of my studies, I’m trying to cluster co-occurrences of URLs and tags in data from Delicious. I found a promising method for this in a paper called “Emergent Semantics from Folksonomies: A Quantitative Study” (pages 6-13). It used a Separable Mixture Model (SMM, described in the paper “Statistical Models for Co-occurrence Data” pages 2-4) to model the data and an adapted EM- algorithm to fit the known data to the model. I coded the algorithm with Java and ran it with a little piece of real data from Delicious. Unfortunately, the results did not seem correct. The results showed that each tag had equal (although varying from tag to tag) possibility to belong to each concept. Now, while this problem could have came from me simply coding the adapted EM- algorithm wrong, I would also like to rule out the possibility of incorrectly initialized variables. This time, since I didn’t know any better way to do it, I simply initialized all the $R_{r\alpha}$ (variables that denote the possibility of co-occurrence $r$ to have raised from concept $\alpha$) to be equal, $1/K$ ($K$ being the number of concepts). My question is two-fold.   Could the flat results come from the flat variable initialization?   How should I initialize the variables from the EM-algorithm in this case?


title: Item set similarity for real-time calculation purposes?

text: I am working on a real-time recommendation engine. At one step, I have a feature that resembles a string-encoded item set, so I took Jaccard on the tokenized string to get a good similarity result. Frankly, Jaccard offers great results, but it takes far too much time on runtime. My current prototype environment for this is written in Java and I am already using the StringTokenizer to fast transform my string into a token set but the set operations thereafter still take a large amount of time - so my idea was to skip the tokenizing part and work on the string directly, but I am not sure on what metric to count for this puspose. Anyone have any ideas or experiences?


title: Logistic regression and clustering?

text: This is regarding a binary classification with logistic regression with one independent variable (IV) and one dependent variable (DV). The IV seems to be distributed normally, the residuals have a random distribution, and AS FAR AS I KNOW, based on my domain knowledge and looking at a plot, it should affect the DV in a linear manner, and it is definitely a monotonic relationship. So basically, it should be a very straight forward relationship. There are 1100 datapoints. When I run the regression and check the predictions with 10fold cross validation, it predicts at about 51% accuracy or even just randomly. There is a lot of variance so the best we can hope for is maybe 54%. It is worth noting the model predicts a success at about 4 to 1 rate when the true rate is close to 1 to 1. **Question 1: At this point, if you were trying to improve this regression, what would you try to do? What do you think may be the problem?** So I tried various things like adding different IV's, trying random transformations, and had no success. Then I started clustering the data and only selected ABOUT half the datapoints which most closely match the test case I am trying to predict. This improved the accuracy greatly. **Question 2: Why would clustering help in such a simple, basic, regression situation? Is there a more philosophically sound, more effective way than clustering to set up this regression?**


title: Questions about Likelihood Principle

text: I currently try to understand Likelihood Principle and I frankly don't get it at all. So, I will write all my question as a list, even if those might be pretty basic questions.   * What exactly does "all of the information" phrase mean in the context of this principle? (as in _all of the information in a sample is contained in the likelihood function.)_   * Is the principle somehow connected to the very provable fact, that $p(x|y)\propto p(y|x)p(x)$? Is the "likelihood" in the principle the same thing, as $p(y|x)$, or not?   * How can a mathematical theorem be "controversial"? My (weak) understanding of math is that a theorem is either proven, or is not proven. To what category does Likelihood Principle fall?   * How is the Likelihood Principle important for Bayesian inference, which is based on $p(x|y)\propto p(y|x)p(x)$ formula?


title: How to decompose a distribution with two peaks?

text: I'm modeling saving account's amount, whose change looks like a log-normal distribution. It means suppose $y$ is total saving account's amount; $x = \ln(y)$ is the natural log; $dx$, the daily change, seems follows normal distribution. However, $dx$ seems has 2 peaks, one is high and a bit less than 0, one is lower and bigger than 0. The reason could be that people are likely to frequently withdraw little amount of money, e.g. weekly withdraw 200 bucks for everyday expense, and deposit more money less frequently, e.g. monthly deposit the paycheck. Now the question is, how could this distribution be decomposed? Seems I can't decompose it to two normal distribution, as two independent normal distribution's sum is still a normal distribution. Any suggestion?


title: GLM model validation; Patterns in residuals ? - R

text: I have attached an image obtained during my model validation. I fitted a negative binomial GLM to my data and I have 4 models possible with more or less the same residual plots. I have been reading about GLM as it is the first time that I use this model and it seems that if there is no pattern in the residuals the model is fine. I have to admit that I cannot really decide if my residuals follow a pattern or not. Is there a test that I could use instead of the graph? If you think that there are patterns here what can I do? I have tried Poisson and quasi poisson as well but there are clear patterns this time. My response variable is a count variable and I have both categorical and continuous covariates. Thank you very much for your help! ![enter image description here](http://i.stack.imgur.com/y0aiq.jpg)


title: Markov random field and iterated condition mode

text: I have spent a lot of time studying MRF (applied to images) but still can't grasp the idea. Could you please clarify these ideas:   1. What is the clique potential? What is a clique in image, and do they overlap?   2. In ICM, people do MAP estimation. As far as I know, posterior is P(parameter| data), so what are parameters here?   3. They also applied Gibbs sampling in estimation. What probability distribution are they sampling? And why you need to do sampling?


title: ANOVA results do not match post-hoc Tukey test, how to proceed?

text: There has been some discussion about why differences exist between factor significance after an ANOVA and no significance in the pairwise comparisons, and it has been very helpful. However, there is no directions regarding how to proceed. For instance, in my case the overall P value given by a one-way ANOVA for the treatment factor was P=0.0127 but all of the treatments came out as not significantly different from each other following Tukey's LSD. Should I run a different post-hoc test in this case? Which one would you recommend? Thank you so much!


title: Question about Dummy Variable in Cross Level Interaction - GENLINEMIXED (SPSS)

text: Hi I am running linear GENLINEMIXED in SPSS 22. What does it mean when one of the dummy variables you are using in a model is significant when one group is coded as a reference group, and when you swap the reference coding to the other group, it is no longer significant. Could it be that one group has a linear relationship to the DV, and the other non-linear? The dummy variable is a level 2 predictor for a level 1 DV. In my case, I have been working with some other co-authors on a paper that is using GENLINMIXED to run a linear multi-level model nested by 54 countries in SPSS on three different dependent variables. These DV variables are ipsative (they add up to 100%). The items ask entrepreneurs and business owners to allocate 100 points across the economic, social, and environmental values of their respective venture to try and capture some quantification of what we call in business research the "triple bottom line." I originally used MANOVA to test the relationships on the DVs simultaneously for this study. However, a reviewer suggested I do a multilevel model for the analysis, and another reviewer I have another reviewer suggesting a Multi-level Seemingly Unrelated Regression. Any suggestions on which route would be best? The following exerpt from our working paper has utilized the multi-level model approach (e.g., GENLINMIXED in SPSS). My goal is to identify a cross-level interaction between gender (coded as binary 0=men 1=women) and culture as operationalized by the World Values Survey (coded binary 0= post-materialist culture and 1=materialist culture). I want to identify if culture attenuates or increases the value goals between the sexes. A weird thing happening with my data when I pursued the GENLINMIXED linear multi-level model is that if I make the reference group for my two binary independent variables men and post-materialist, the cross level interaction for women and materialism is significant across all three DVs (economic, social, and environmental). Also, the main effect of gender is significant, and culture not significant. However, if I make the reference group for the binary independent variables women and materialism, the cross level interaction term men and post-materialisim is not significant. Again gender is significant and culture not significant. Again, could it be that the cross level interaction is linearly related to the DVs for women in materialist countries, and not linearly related to the DVs for men in post-materialist countries? I have just the methods write-up of that paper attached two with the three hypothesis if you would like further detail into the situation. Any advice or feedback would be most appreciated.


title: Pros and cons of meta-analyses

text: I have been considering doing some meta-analysis for a particular field of study in evolution, but before I go any further I would like to know; what are the positives and negatives of the process? For example, no need for a practical experiment is an advantage (time & money) but there will be a publication bias (more exciting results get published) which would be a disadvantage. What papers in statistics journals discuss the pros and cons of meta-analysis?


title: Permutation test for the significance of several dependent regression coefficients from separate models

text: I have $i>100$ Cox regression models which differ only in one explanatory variable, $x_i$. They have the same outcome data and the same other explanatory (adjustment) variables. I want to control the global Type 1 error probability/adjust for multiple testing of the significance of the coefficients of the $x_i$. Some of the $x_i$ show high pairwise correlations. I inferred that therefore their coefficients and p-values are dependent. Neither Bonferroni nor Holm correction take into account these dependencies. The most attractive/exact option in my eyes is a permutation test that does consider them. I thought about scrambling the rows of the outcome survival matrix, thereby conserving the correlations between the $x_i$ while destroying their relationship with the outcome. Then I compute the coefficients of all the $x_i$. In doing this for a few thousand different resamples I get the multivariate null distribution of the coefficients. And this is where I'm stuck. I wanted to take the proportion of the null/resampled distribution for each coefficient estimated in the original data that is higher than it (or lower in the case of negative coefficients) as its p-value. But that would, again, destroy the dependency structure, as it treats the multivariate distribution as $i$ separate univariate ones. How would one go about doing a permutation test for this problem? Is it even possible?


title: practical vs statistical significance

text: In a nutshell what would be the main differences between practical and statistical significance? What would be the relative importance of each in drawing scientific inferences, for example?


title: How to calculate the time 2 standard deviation given time 1 and change standard deviations?

text: I need to calculate the mean(+/-SD) final value of a variable from the mean(+/-SD) baseline value and the mean(+/-SD) change from baseline. So, for example:   * Mean baseline weight is 40(+/-2) kg.   * Mean change from baseline is +5(+/-1) kg. I'd assume the correct final mean is 45kg, but how would I go about calculating a new SD? I've googled for methods of calculating 'Pooled SD', but those only seem to apply to pooling same variables (like weights from different groups) and I don't think those apply here. This question seems the most similar to my question but also seems to try the exact opposite of what I'm trying to achieve.


title: Issues with tritomising and dichotomising variable

text: I see that the issue of splitting a variable into two categories has been discussed on this forum and I understand pros, cons as well as some ways of performing it. The issue that I am interested and can't seem to find anything on splitting a variable into three points. For this to make sense let me briefly explain my data. I have used a psychographic scale (5 Point Likert) which consists of eight items where low scores indicate one trait and high scores indicate a different one. Theoretically, most people should fall into the mid group, with much smaller extreme groups which show clear tendencies for certain behaviour. Thus, I used this scale with an intention of grouping my respondents into low, middle, and high groups. This classification into three groups is supported conceptually by previous research but no study discloses how groups are split. I tried using a composite score for the scale and e.g. correlate with a dependent variable. However, because extreme groups are relatively small the result for the whole sample indicates asssociation in the expected direction albeit very low. I have considered dichotomising but individuals with top or low scores on the eight variables get muddled up with middle scores that conceptually do not indicate any distinct group. Having looked into graphical representation of the scale I can conclude that distribution is normal and indeed most people score in the middle. Thus, visually I can identify the cut-off points which seems to make sense, it produces group sizes supported theoretically and any further tests show significant differences between groups on dependent variables of interest. However, I understand that cutting data this way is very arbitrary. Is there any way I can actually justify doing this based on conceptual support for groups? Is there anything else anyone can suggest as I have run out of ideas. Any help will be much appreciated. Many thanks in advance.


title: Linear model question regarding scale of variables

text: I want to perform linear regression using the following command: `lm(clinicVisits ~ clinicRank)` where `clinicVisits` is the number of visits to a clinic, and `clinicRank` is the rank of the clinic. `clinicVisits` is an integer in the range of 1–200K, and clinicRank is a continuous variable and has a range of 1.0–150.0 where 1.0 is the best ranking. I am trying to predict what an increase in rank will mean in terms of number of visits, as I believe there is a relationship between the two, so I am seeking this linear regression equation:               clinicVisits = y intercept + clinicRank*(b1)      I also assume that `b1` will be an increasingly negative number as the `clinicRank` gets smaller (i.e. closer to 1, the best rank). Do I need to apply a scale to the variables since the `clinicRank` is relatively small compared to visits, or a log to the `clinicVisits` since it has a relative large range?


title: Choice of 0 or -1 for failure in the independent variables of a logistic regression

text: I am performing some exploratory analysis on a dataset where the dependent variable is a dichotomous variable. I have ~10 explanatory variables, some of which are dichotomous observations. I am running a logistic regression to learn more about the dataset. When converting observations to quantitative measures I always select 1 as the success state (e.g. the user successfully completed a funnel). I am looking to learn whether I should choose 0 or -1 for the failure state. I am not certain of the statistical implications of the choice - what are they? How would the choice affect the analysis? Does the choice of 0 not conflate with missing data?


title: Can we approximate this function by a polynomial?

text: For the more mathematically minded, we have $x \in \mathbb{R}^2$ and the function $h(x)$ defined as: $h(x)=\alpha_1x_1^2+\alpha_2x_2^2+\alpha_3x_1+\alpha_4x_2+\alpha_5x_1x_2+\alpha_6$ and the vector of alpha's is **known** and further guaranteed to be such that $h(x)$ is a general elliptic paraboloid (i.e. a convex function). Now, define $g(x)=(h(x))^+$ where $(z)^+$ is the positive part of $z\in\mathbb{R}$. The questions follows:   1. Is there a way to globally approximate $g(x)$ by a polynomial ? If so what it is?   2. Can this approach (the answer to the previous sub-question) be extended to $x\in\mathbb{R}^p$ with $p$ moderatly large ?   3. Would the problem be any easier if we were to assume $\alpha_5=0$ ? Following Whuber's comment: is it possible to find a polynomial approximation to $g(x)$ that would be better than $h(x)$? in the event that many approximations solution for this problem exist, what are they ? I can obviously solve $argmin.\int_{\mathbb{R}^2}(g(x)-\hat{g}(x))^2 dx$. I'm wondering if there are explicit, known solutions to this problem (i.e. polynomial series for example of which i know close to nothing).


title: Which variables should I transform, center, and/or standardize in my data for Principal Component Analysis?

text: I have multiple datasets that I am attempting to use principal component analysis (PCA) on in order to infer the underlying structure of the data. I'm attempting to predict growth increment of conifers (the response which apparently needed to be cube root transformed for my original linear models). My predictor variables include: continuous (some log-transformed like size- related variables, some compositional, and some traditionally scaled); binary; ordinal; and nominal (of which I'm pretty sure I can order with respect to site productivity). I'm mostly attempting to make inferences because some of my previous models from linear regression have not produced significant terms corresponding to categories like health or climate and I feel that they should have. I've read in Joliffe (2002) that I might want to leave my compositional data uncentered. Forgive my naiveté, but I'm wondering if I should also do this for my binary data. It also seems that I might not want to scale these since they're already on a 0 - 1 scale but I'm not sure. After only using my continuous predictors (including compositional variables), it appears that a size principal component makes up the largest portion of the variability and Joliffe also mentions that one might want to attempt a "double-centering" of the log-transformed data if this is the case. I don't know if this is even possible since the log-transformed data only makes up a small portion of my datasets. Finally, I should ask if having such a combination of data types will complicate analysis so much that it will make inferences unreliable. I am currently using R for this analysis. Thank you.


title: Pitman's test of equality of variance and testing for regression to the mean: am I doing the right thing?

text: I have 2560 paired observations from an experiment in which participants provided two ratings for a set of objects, at two different points in time. Half of the objects in the set had the value of an attribute A changed between the two time points, half did not. Of the objects that were changed in each participant's set, half went from A' to A'' and half from A'' to A'. (i.e. all participants experienced both orders). My main hypothesis is that changing this attribute from A' to A'' leads on average to a higher rating, and this is indeed supported by the data. I am also interested in determining whether the magnitude and perhaps direction of this effect depends on the A' rating. For the purposes of this question, I am considering only those instances where A was changed (1280 pairs of obs). The following GLMM (A'' rating - A' rating) = participant + order + A' rating where A' rating is a covariate and participant and order are categorical variables, leads to the conclusion that there is a significant positive correlation between A' rating and effect of changing to A'' and that this correlation is <1, such that objects with a low A' rating have their rating increased by changing to A'' but that objects with a high A' rating actually get rated lower when changed to A''. I want to test whether this is simply due to regression to the mean. To this end, I have followed Kelly and Price in using Pitman's test of equality of variances for paired samples and would appreciate some feedback on whether I've done the right thing. This is what I did, following the suggestion of a colleague: 1) calculated the SD of A'' ratings $(SD_1)$ and the SD of A' ratings $(SD_2)$   2) regressed A'' rating on A' rating and recorded the correlation $r$.   3) calculated T as $T=\frac{\sqrt{(n-2)} [(SD_1/SD_2)-(SD_2/SD_1)]}{2 \sqrt{(1-r^2)}}$ The 2-tailed p value of T (Student's t dist with 1280-2 DF) is 0.07, i.e. at alpha=0.05 there is no significant difference between the variances for the two sets of ratings and thus no effect of A' rating on rating difference beyond regression to the mean. (We can argue about 2-t vs. 1-t p values later). I now plan to adjust my difference scores to account for this and re-do the GLMM outlined above, as outlined by Kelly & Price. If you've got this far through the detail, then firstly well done, and secondly, can you tell me if 1) my colleague's suggestion was sensible and 2) if it was, have I implemented it correctly? I have a couple of concerns/apparent grey areas but I'd be interested to hear what others have to say first. Thank you.


title: Bayes Estimation on Dirichlet distribution

text: I'm trying to get my head around the "hidden species" problem. It goes something like this. You visit a park and run into three species, 3 lions, 2 tigers and 1 bear. You are to determine what is your best guess on the distribution of these species. You may be tempted to do $\frac{3}{3+2+1} = 50\%$ for the lions and 33% for the tigers 17% for the bears, but this is wrong. Instead the dirichlet distribution is used to solve this. Described here. In short, what the author says is to use $\alpha_{i} = \\{1,1,1\\}$ and arrive at $P(Lion) = \frac{1 + 3}{1+1+1 + 3 + 2 + 1} = \frac{4}{9}$. The $\\{1,1,1\\}$ in the denominator is the sum of the $\alpha_i's$ Specifically my questions are the following. 1) While I get the fact that the dirichlet distribution is the conjugate prior multinomial distribution, how did the author choose the alpha's in the dirichlet distribution to be (1,1,1)? 2) Since the choice of alpha's here are like "imaginary" counts, how can I accomodate the info (if it exists) that I'd seen 4 lions, 1 tiger & 2 bears (say) in the past. Could someone shed some more explanatory light on this? Thanks much